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ABSTRACT 

Fermi  has  provided  the  largest  sample  of  y -ray- selected  blazars  to  date.  In  this  work  we  use  a  complete  sample  of 
flat  spectrum  radio  quasars  (FSRQs)  detected  during  the  first  year  of  operation  to  determine  the  luminosity  function 
(LF)  and  its  evolution  with  cosmic  time.  The  number  density  of  FSRQs  grows  dramatically  up  to  redshift  ~0. 5-2.0 
and  declines  thereafter.  The  redshift  of  the  peak  in  the  density  is  luminosity  dependent,  with  more  luminous  sources 
peaking  at  earlier  times;  thus  the  LF  of  y- ray  FSRQs  follows  a  luminosity-dependent  density  evolution  similar 
to  that  of  radio-quiet  active  galactic  nuclei.  Also,  using  data  from  the  Swift  Burst  Alert  Telescope  we  derive  the 
average  spectral  energy  distribution  (SED)  of  FSRQs  in  the  10  keV-300  GeV  band  and  show  that  there  is  no 
correlation  between  the  luminosity  at  the  peak  of  the  y-ray  emission  component  and  its  peak  frequency.  Using 
this  luminosity-independent  SED  with  the  derived  LF  allows  us  to  predict  that  the  contribution  of  FSRQs  to  the 
Fermi  isotropic  y- ray  background  is  9.3+}{60%  (±3%  systematic  uncertainty)  in  the  0.1-100  GeV  band.  Finally  we 
determine  the  LF  of  unbeamed  FSRQs,  finding  that  FSRQs  have  an  average  Lorentz  factor  of  y  —  1L7*232,  that 
most  are  seen  within  5°  of  the  jet  axis,  and  that  they  represent  only  ~0.1%  of  the  parent  population. 

Key  words:  cosmology:  observations  -  diffuse  radiation  -  galaxies:  active  -  galaxies:  jets  -  gamma  rays:  diffuse 
background  -  surveys 
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1.  INTRODUCTION 

The  detection  of  luminous  quasars  at  redshift  >6  (e.g., 
Fan  et  al.  2003;  Willott  et  al.  2010)  provides  evidence  of 
supermassive  black  hole  (SMBH)  formation  in  the  first  1  Gyr  of 
cosmic  time.  There  are  appreciable  challenges  to  forming  (see, 
e.g.,  Wyithe  &  Loeb  2003;  Volonteri  &  Rees  2005;  Begelman 
et  al.  2006;  Mayer  et  al.  2010)  and  fueling  (see  Kauffmann 
&  Haehnelt  2000;  Wyithe  &  Loeb  2003;  Croton  et  al.  2006) 
these  objects  at  such  early  times,  although  it  is  widely  believed 
that  strong  accretion  can  be  initiated  by  major  mergers  (see 
Kauffmann  &  Haehnelt  2000;  Wyithe  &  Loeb  2003;  Croton 
et  al.  2006). 

Blazars  represent  an  extreme  manifestation  of  such  active 
galactic  nucleus  (AGN)  activity,  with  radiation  along  Earth’s  line 
of  sight  dominated  by  a  relativistic  jet.  It  is,  as  yet,  unclear  how 
such  jet  activity  connects  with  the  more  isotropically  emitted 
bulk  accretion  luminosity.  For  example  according  to  Blandford 
&  Znajek  (1977),  the  energy  stored  in  a  black  hole’s  spin  can  be 
extracted  in  the  form  of  a  relativistic  jet.  Thus,  blazar  evolution 
may  be  connected  with  the  cosmic  evolution  of  the  spin  states 
of  massive  black  holes.  Radio-loud  (RL,  jet  dominated)  blazars 
have  been  seen  at  redshifts  as  high  as  z  =  5.5  (Romani 
2006),  and  it  is  plausible  that  major  mergers,  more  frequently 
experienced  in  the  early  universe,  might  preferentially  produce 
maximally  rotating  black  holes  (e.g.,  see  Escala  et  al.  2004; 
Dotti  et  al.  2007). 

Thus  the  study  of  RL  AGN,  blazars,  with  strong  relativisti- 
cally  beamed  jets  can  provide  a  method  to  study  jet  activity, 
BH  spin,  and  major  merger  events.  This  can  be  done  by  de¬ 
termining  the  luminosity  function  (LF)  of  blazars  (essentially 
the  number  of  blazars  per  comoving  volume  element  within  a 


certain  luminosity  range)  and  its  evolution  with  redshift.  The 
Fermi  Gamma-ray  Space  Telescope  provides  one  of  the  largest 
data  sets  with  which  to  study  the  properties  of  blazars.  Thanks 
to  its  sensitivity  and  uniform  coverage  of  the  sky,  Fermi  has 
detected  hundreds  of  blazars  from  low  redshifts  out  to  z  =  3.1 
(Ackermann  et  al.  201 1). 

The  LF  of  blazars  also  allows  us  to  evaluate  their  contribution 
to  the  diffuse  backgrounds  and  to  determine  their  relationship 
with  the  parent  population  (Ajello  et  al.  2008b;  Inoue  2011). 
Blazars  have  been  extensively  studied  at  radio  (Dunlop  & 
Peacock  1990;  Wall  etal.  2005),  soft  X-ray  (Giommi  &  Padovani 
1994;  Rector  et  al.  2000;  Wolter  &  Celotti  2001;  Caccianiga 
et  al.  2002;  Beckmann  et  al.  2003;  Padovani  et  al.  2007),  and 
GeV  energies  (Hartman  et  al.  1999).  It  seems  clear  that  flat 
spectrum  radio  quasars  (FSRQs)  evolve  positively  (i.e.,  there 
were  more  blazars  in  the  past,  Dunlop  &  Peacock  1990)  up  to 
a  redshift  cutoff  which  depends  on  luminosity  (e.g.,  Padovani 
et  al.  2007;  Wall  2008;  Ajello  et  al.  2009a).  In  this  respect 
FSRQs  evolve  similarly  to  the  population  of  X-ray-selected, 
radio-quiet  AGNs  (Ueda  et  al.  2003;  Hasinger  et  al.  2005;  La 
Franca  et  al.  2005).  On  the  other  hand,  the  evolution  of  the  other 
major  class  of  Fermi- detected  AGN,  BL  Lac  objects,  and  their 
relation  to  FSRQs,  remains  a  matter  of  debate,  with  claims  of  no 
evolution  (Caccianiga  et  al.  2002;  Padovani  et  al.  2007)  or  even 
negative  evolution  (e.g.,  Rector  et  al.  2000;  Beckmann  et  al. 
2003).  Samples  with  larger  redshift  completeness  fractions  are 
needed  to  validate  these  claims. 

In  this  work,  we  report  on  the  LF  of  FSRQs  detected  by 
Fermi  in  its  first  year  of  operation.  There  have  been  attempts  in 
the  past  (e.g.,  Chiang  et  al.  1995)  to  characterize  the  evolution  of 
y- ray  AGNs,  starting  from  the  EGRET  sample  (Hartman  et  al. 
1999).  One  challenge  was  the  small  sample  size  and  redshift 
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incompleteness  of  the  EGRET  set.  Often  it  was  assumed  that 
the  y  -ray-detected  blazars  had  an  LF  following  that  of  another 
band,  e.g.,  radio-  or  X-ray- selected  blazars.  The  results  reported 
in,  e.g.,  Stecker  &  Salamon  (1996),  Narumoto  &  Totani  (2006), 
Inoue  &  Totani  (2009),  and  Stecker  &  Venters  (2011)  follow 
this  approach.  Alternatively,  an  LF  may  be  estimated  from  the 
]/-ray  sample  directly  (Chiang  &  Mukherjee  1998;  Miicke  & 
Pohl  2000;  Dermer  2007;  Bhattacharya  et  al.  2009),  although 
only  a  small  (~60)  blazar  sample  including  both  BL  Lac  objects 
and  FSRQs  was  available  (from  EGRET  data)  with  acceptable 
incompleteness. 

We  report  here  on  a  detailed  LF  measured  from  a  sample 
of  186  y  -ray- selected  FSRQs  detected  by  Fermi.  The  work  is 
organized  as  follows.  Sections  2  and  3  describe  the  properties  of 
the  sample  used  and  the  method  employed  to  determine  the  LF 
of  blazars.  The  LF  of  FSRQs  is  derived  in  Section  4.  In  Section  5, 
the  spectral  energy  distributions  (SEDs)  of  Fermi's  FSRQs  are 
analyzed  in  detail,  testing  for  possible  correlations  of  the  peak 
energy  with  the  peak  luminosity.  The  contribution  of  FSRQs  to 
the  isotropic  gamma-ray  background5  (IGRB;  see  Abdo  et  al. 
2010e)  is  determined  and  discussed  in  Section  6.  Throughout 
this  paper,  we  assume  a  standard  concordance  cosmology  (Hq  = 
71  km  s-1  Mpc-1  and  QM  =  1—  Qa  =  0.27). 

2.  THE  SAMPLE 

The  First  Fermi  Large  Area  Telescope  (LAT)  Catalog  (1FGL; 
Abdo  et  al.  2010a)  reports  on  more  than  1400  sources  detected 
by  Fermi -LAT  during  its  first  year  of  operation.  The  first  LAT 
AGN  catalog  (1LAC;  Abdo  et  al.  2010g)  associates  ~700  of  the 
high-latitude  1FGL  sources  (\b\  ^  10°)  with  AGNs  of  various 
types,  most  of  which  are  blazars.  The  association  method  is 
probabilistic,  by  comparing  the  chance  probability  for  gamma- 
ray  sources  to  be  found  in  the  vicinity  of  sources  in  some 
underlying  radio,  optical,  and  X-ray  AGN  catalogs  with  the 
probability  of  chance  associations.  Among  these  catalogs,  the 
most  important  ones  are  the  Combined  Radio  All-sky  Targeted 
Eight  GHz  Survey  (CRATES;  Healey  et  al.  2007),  the  Candidate 
Gamma-Ray  Blazar  Survey  (CGRaBS;  Healey  et  al.  2008),  and 
the  Roma-BZCAT  (Massaro  et  al.  2009),  which  is  a  catalog  of 
blazars  identified  at  all  frequencies.  Here,  we  adopt  the  1LAC 
associations  with  probability  >0.8,  i.e.,  <20%  chance  of  being 
a  false  positive,  although  most  of  the  identifications  are  at  the 
0.95  or  better  level  (Abdo  et  al.  2010g).  The  subset  of  these 
1LAC  AGNs  used  in  the  present  analysis  are  those  sources  at 
\b\  >  15°  detected  by  the  pipeline  developed  by  Abdo  et  al. 
(201  Of)  with  a  test  statistic  (TS)  significance  greater  than  50 
(corresponding  to  a  significance  of  >7a;  see,  e.g.,  Mattox  et  al. 
1996).  For  this  sample,  we  have  produced  a  set  of  Monte  Carlo 
simulations  that  can  be  used  to  determine  and  account  for  the 
selection  effects.  This  sample  contains  483  objects,  1 86  of  which 
are  classified  as  FSRQs.  The  faintest  identified  FSRQ  has  a 
100  MeV-100  GeV  band  flux  of  Fm  ~  10-8  photons  cm-2  s-1 . 
To  limit  the  incompleteness  (i.e.,  the  fraction  of  sources  without 
an  association)  we  apply  this  as  a  flux  limit  to  the  sample 
of  483  objects,  resulting  in  a  full  sample  of  433  sources  of 
which  29  (i.e.,  ~7%)  do  not  have  associated  counterparts.  Note 
that  all  of  the  associations  with  an  FSRQ  designation  have 
spectroscopically  measured  redshifts  either  from  the  literature 
or  from  recent  spectroscopic  studies  (Shaw  et  al.  2012). 

5  The  isotropic  gamma-ray  background  refers  to  the  isotropic  component  of 

the  Fermi  sky  (Abdo  et  al.  2010e)  and  as  such  might  include  components 

generated  locally  (e.g.,  Keshet  et  al.  2004)  and  components  of  truly 

extragalactic  origin. 


Tabic  1 

Composition  of  the  \b\  ^  15°,  TS  ^  50,  Fioo  ^  10-8  photons  cm-2  s-1 
Sample  Used  in  This  Analysis 


CLASS 

No.  of  objects 

Total 

433 

FSRQs 

186 

BL  Lac  objects 

157 

Pulsars 

28 

Othera 

16 

Radio  associations'3 

17 

Unassociated  sources 

29 

Notes. 

a  Includes  starburst  galaxies,  LINERS,  narrow-line  Seyfert  1  objects, 
and  Seyfert  galaxy  candidates. 

b  Fermi  sources  with  a  radio  counterpart,  but  no  optical  type  or 
redshift  measurement. 

The  composition  of  this  sample  is  reported  in  Table  1.  The 
186  FSRQs  detected  by  Fermi  with  TS  ^  50,  \b\  >  15°,  and 
Tqoo  ^  10-8  photons  cm-2  s-1  constitute  the  sample  that  will 
be  used  in  this  analysis. 

3.  METHOD 

A  classical  approach  to  derive  the  LF  is  based  on  the 
1/Vmax  method  of  Schmidt  (1968)  applied  to  redshift  bins. 
However,  this  method  is  known  to  introduce  bias  if  there  is 
significant  evolution  within  the  bins.  Moreover,  given  our 
relatively  small  sample  size  and  the  large  volume  and  luminosity 
range  spanned,  binning  would  result  in  a  loss  of  information. 
Thus,  we  decided  to  apply  the  maximum-likelihood  (ML) 
algorithm  first  introduced  by  Marshall  et  al.  (1983)  and  used 
recently  by  Wall  et  al.  (2008)  and  Ajello  et  al.  (2009a)  for 
the  study  of  (respectively)  submillimeter  galaxies  and  blazars 
detected  by  Swift.  The  aim  of  this  analysis  is  to  determine  the 
space  density  of  FSRQs  as  a  function  of  rest-frame  0. 1-100  GeV 
luminosity  (Ly),  redshift  (z),  and  photon  index  (F),  by  fitting  to 
the  functional  form 

d3N  d2N  dN  dV  dN  dV 

dLydzdT  dLydV  dr  dz  7  dr  dz 

(1) 

where  <D(LK ,  z)  is  the  LF,  and  dV /dz  is  the  co-moving  volume 
element  per  unit  redshift  and  unit  solid  angle  (see,  e.g.,  Hogg 
1999).  The  function  dN/dT  is  the  (intrinsic)  photon  index 
distribution  and  is  assumed  to  be  independent  of  z.  It  is  modeled 
as  a  Gaussian: 

dN  (x-n)2 

—  (2) 
dr 

where  ji  and  a  are,  respectively,  the  mean  and  the  dispersion 
of  the  Gaussian  distribution.  The  spectrum  of  most  blazars  is 
known  to  deviate  from  a  simple  power  law  (see,  e.g.,  Abdo 
et  al.  2010d,  and  Section  5).  However,  detection  of  sources  in 
Fermi  is  performed  by  means  of  an  ML  fit  (to  the  source  and 
background  photons)  by  modeling  the  source  spectrum  with  a 
power  law  (Abdo  et  al.  2010a).  For  consistency,  we  follow  here 
the  same  procedure  and  assess  in  the  Appendix  the  uncertainties 
connected  with  this  hypothesis. 

The  best-fit  LF  is  found  by  comparing,  through  an  ML 
estimator,  the  number  of  expected  objects  (for  a  given  model  LF) 
to  the  observed  number  while  accounting  for  selection  effects  in 
the  survey.  In  this  method,  the  space  of  luminosity,  redshift,  and 
photon  index  is  divided  into  small  intervals  of  size  dLydzdr.  In 
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each  element,  the  expected  number  of  blazars  with  luminosity 
Ly ,  redshift  z,  and  photon  index  T  is 

k(Ly,z,  T)dLydzdT  =  O (Ly,  z)Q(Ly,  z,  T) 

dN  dV 

x  ————dLydzdr,  (3) 
dT  dz 

where  £2(LK,z,r)  is  the  sky  coverage  and  represents  the 
probability  of  detecting  in  this  survey  a  blazar  with  luminosity 
Ly ,  redshift  z,  and  photon  index  T.  This  probability  was  derived 
for  the  sample  used  here  by  Abdo  et  al.  (201  Of)  and  the  reader 
is  referred  to  that  aforementioned  paper  for  more  details.  With 
sufficiently  fine  sampling  of  the  Ly  —  z— T  space  the  infinitesimal 
element  will  either  contain  0  or  1  FSRQs.  In  this  regime  one  has 
a  likelihood  function  based  on  joint  Poisson  probabilities: 

L  =  tl  MLy  i,  zi,  ri)dLYdzdre-'ML''-i’z‘Ji>dLYdzdr 

i 

X  J- J  g—^(Lyj,Zj,Tj)dLrdzdr 

j 


This  is  the  combined  probability  of  observing  one  blazar  in 
each  bin  of  (Lyj,  Zi,  F/)  populated  by  one  Fermi  FSRQ  and 
zero  FSRQs  for  all  other  (Lyj,  zj,  Fj).  Transforming  to  the 
standard  expression  B  =  —  2  In  L  and  dropping  terms  which 
are  not  model  dependent,  we  obtain 


B  =  -  2  In 


d3N 


dLydzdT 


+  2 


P  T'rnax  p^y, max  pZ 

7  Fnin  Jl'v.min  7  Zm 


X(Ly ,  T,  z)dLydzdr.  (5) 


dN 

dLy 


A .(Ly,  T,  z)dzdr 


(7) 


dN 

~dr 


A .(Ly ,  T,  z)dLYdz , 


(8) 


where  the  extremes  of  integrations  are  the  same  as  in 
Equation  (5).  Since  Ly  depends  on  redshift,  Equations  (6) 
and  (7)  are  not  independent.  The  source  count  distribution  can 
be  derived  as 


N(>S)  = 


P  Fnax  pZmax  p^y,n 
7  Train  J  Zmin  7  Ly{z, 


dN  dV 

<r>(Ly,z)-——-drdzdLy, 
S)  dr  dz 


(9) 


where  Ly(z,  S)  is  the  luminosity  of  a  source  at  redshift  z  having 
a  flux  of  S. 

To  display  the  LF  we  rely  on  the  ‘Wobs/ Amdl”  method  devised 
by  La  Franca  &  Cristiani  (1997)  and  Miyaji  et  al.  (2001)  and 
employed  in  several  recent  works  (e.g.,  La  Franca  et  al.  2005; 
Hasinger  et  al.  2005).  Once  a  best-fit  function  for  the  LF  has 
been  found,  it  is  possible  to  determine  the  value  of  the  observed 
LF  in  a  given  bin  of  luminosity  and  redshift: 


Zi)  =  Omdl(Lri,-,  Zi) 


N°hs 
Ymdl  ’ 


(10) 


where  Lyj  and  Zi  are  the  luminosity  and  redshift  of  the  ith  bin, 
O mdl(Lyj ,  Zi )  is  the  best-fit  LF  model,  and  A°bs  and  Nfdl  are  the 
observed  and  the  predicted  number  of  FSRQs  in  that  bin.  These 
two  techniques  (the  Marshall  et  al.  1983  ML  method  and  the 
“7Vobs/7Vmdl”  estimator)  provide  a  minimally  biased  estimate  of 
the  LF  (cf.  Miyaji  et  al.  2001). 


The  limits  of  integration  of  Equation  (5),  unless  otherwise  stated, 
are  Ly  m in  =  10  erg  s  ,  Lymax  =  106  erg  s  ,  Zmin  = 
10-2,  Zmax  =  6,  rmin  =  1.8,  and  rmax  =  3.0.  The  results  of 
the  analysis  do  not  depend  on  the  limits  of  integration  zmax 
and  Ly,max.  The  values  of  the  Ly^m jn  and  zmin  are  chosen  to 
be  a  factor  of  a  few  lower  than  the  smallest  values  observed 
in  Ackermann  et  al.  (2011)  to  force  the  LF  to  account  for 
the  paucity  of  low-luminosity  low-redshift  FSRQs.  However, 
we  get  results  compatible  within  the  statistical  uncertainties 
if  we  use  the  minimum  observed  luminosity  and  redshift  of 
the  source  sample  of  Table  1.  The  best-fit  parameters  are 
determined  by  minimizing6  B  and  the  associated  la  error 
is  computed  by  varying  the  parameter  of  interest,  while  the 
others  are  allowed  to  float,  until  an  increment  of  A B  =  1  is 
achieved.  This  gives  an  estimate  of  the  68%  confidence  region 
for  the  parameter  of  interest  (Avni  1976).  While  computationally 
intensive,  Equation  (5)  has  the  advantage  that  each  source  has 
its  appropriate  individual  detection  efficiency  and  ^-correction 
treated  independently. 

In  order  to  test  whether  the  best-fit  LF  provides  a  good 
description  of  the  data  we  compare  the  observed  redshift, 
luminosity,  index,  and  source  count  distributions  against  the 
prediction  of  the  LF.  The  first  three  distributions  can  be  obtained 
from  the  LF  as 


4.  RESULTS 

4.1.  Pure  Luminosity  Evolution  and  the  Evidence 
for  a  Redshift  Peak 

The  space  density  of  radio-quiet  AGNs  is  known  to  be 
maximal  at  intermediate  redshift.  The  epoch  of  this  “redshift 
peak”  correlates  with  source  luminosity  (e.g.,  Ueda  et  al.  2003; 
Hasinger  et  al.  2005).  This  peak  may  represent  the  combined 
effect  of  SMBH  growth  over  cosmic  time  and  a  falloff  in  fueling 
activity  as  the  rate  of  major  mergers  decreases  at  late  times.  To 
test  whether  such  behavior  is  also  typical  of  the  LAT  FSRQ 
population,  we  perform  a  fit  to  the  data  using  a  pure  luminosity 
evolution  (PLE)  model  of  the  form 

&(Ly,z)  =  0(Ly/e(z)),  (ID 

where 

0(LY/e(z  =  0))  =  ^ 

dLy 

L 


ln(10)Lv 


(hL 

L* 


dN 

dz 


n: 

”  A  min  ”  E  v 


l(Ly,  r,  z)dLydT 


(6) 


6  The  MINUIT  minimization  package,  embedded  in  ROOT  (root.cern.ch), 
has  been  used  for  this  purpose. 


and 

e(z)  =  (1  +z)V/?.  (13) 

In  this  model,  the  evolution  is  entirely  in  luminosity,  i.e., 
the  FSRQ  were  more  luminous  in  the  past  if  positive  evolution 
( k  >  0)  is  found  (the  opposite  is  true  otherwise).  It  is  also 
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Figure  1.  Redshift  (upper  left),  luminosity  (upper  right),  photon  index  (lower  left),  and  source  count  (lower  right)  distributions  of  LAT  FSRQs.  The  line  shows  the 
best-fit  PLE  model  convolved  with  the  selection  effects  of  Fermi  discussed  in  the  text. 


Table  2 

Best-fit  Parameters  of  the  Pure  Luminosity  Evolution  LF 


Sample 

No.  of  Objects 

Aa 

Yl 

L* 

Y2 

k 

5 

U 

Cf 

ALL 

186 

5.59(±0.41)  x  103 

0.29  ±  0.53 

0.026  ±  0.066 

1.25  ±0.32 

5.70  ±  1.02 

-0.46  ±  0.01 

2.45  ±0.13 

0.18  ±0.01 

Low  L 

89 

15.4(±0.2)  x  103 

0.29 

0.026 

1.25 

4.30  ±  2.39 

-0.50  ±  0.04 

2.47  ±  0.04 

0.21  ±  0.03 

High  L 

97 

22.6(±2.0)  x  103 

0.29 

0.026 

1.25 

3.47  ±  1.73 

-0.79  ±  0.04 

2.46  ±  0.02 

0.20  ±  0.02 

Notes.  Parameters  without  an  error  estimate  were  kept  fixed  during  the  fit. 
a  In  units  of  10“ 13  Mpc-3  erg-1  s. 

straightforward  to  demonstrate  that  the  luminosity  evolution 
(i.e.,  Equation  (13))  of  FSRQs  peaks  at  zc  =  —  1  —  k%.  The  best- 
fit  parameters  are  reported  in  Table  2.  The  evolution  of  the  FSRQ 
class  is  found  to  be  positive  and  fast  ( k  =  5.70  zb  1.02).  The 
redshift  peak  is  zc  =  1 .62  zb  0.03.  Moreover,  the  subsequent  rate 
of  decrease  of  the  luminosity  after  the  peak  is  well  constrained 
(§  =  —0.46  zb  0.01).  However,  as  shown  in  Figure  1,  while 
this  model  provides  a  good  fit  to  the  observed  redshift  and 
luminosity  distributions,  it  is  a  very  poor  representation  of  the 
measured  logN-logS. 

We  next  test  whether  the  redshift  peak  depends  on  luminosity, 
splitting  the  data  set  at  Ly  =  3.2  x  1047  erg  s_1.  A  fit  is  then 
performed  to  each  sub-sample  to  determine  the  position  of  the 
redshift  peak  (if  any),  keeping  the  parameters  of  Equation  (12) 
fixed.  The  results  of  the  fits  to  the  low-  and  high-luminosity 
data  sets  are  reported  in  Table  2.  From  Equation  (13)  and  the 
values  of  k  and  §  it  is  apparent  that  there  is  a  significant  shift 
in  the  redshift  peak,  with  the  low-  and  high-luminosity  samples 
peaking  at  ~  1.1 5  and  ~1.77,  respectively. 


4.2.  The  Luminosity -dependent  Density  Evolution 
and  the  Redshift  Peak 

Since  a  simple  PEE  LF  model  provides  an  inadequate  fit  to 
the  Fermi  data  and  since  there  is  some  evidence  for  the  evolution 
of  the  redshift  peak  with  luminosity,  we  now  fit  the  Fermi  FSRQ 
set  to  a  luminosity-dependent  density  evolution  (LDDE)  model. 
Here,  the  evolution  is  primarily  in  density  with  a  luminosity- 
dependent  redshift  peak.  The  LDDE  model  is  parameterized  as 


d>(Ly,z)  =  <I>(Ly)xe(z,Lyf 

where 


e(z,  Ly) 


i+z  y1  /  i+z  y2' 

1  +  Zc(Ly)J  +\l+Zc(Ly)j 


and 

Zc(Ly)  =  Z*  ■  (Ly/1048)“. 


(14) 


(15) 


(16) 
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Figure  2.  Redshift  (upper  left),  luminosity  (upper  right),  photon  index  (lower  left),  and  source  count  (lower  right)  distributions  of  LAT  FSRQs.  The  dashed  line  is  the 
best-fit  LDDE  model  convolved  with  the  selection  effects  of  Fermi.  Note  the  improved  source  count  distribution  over  the  predictions  of  the  PLE  model  of  Figure  1. 


O (Ly)  is  the  same  double  power  law  used  in  Equation  (12). 
This  parameterization  is  similar  to  that  proposed  by  Ueda  et  al. 
(2003),  but  is  continuous  around  the  redshift  peak  zc(Ly). 
This  has  obvious  advantages  for  fitting  algorithms  that  rely 
on  the  derivatives  of  the  fitting  function  to  find  the  minimum. 
Here  zc(Ly)  corresponds  to  the  (luminosity-dependent)  redshift 
where  the  evolution  changes  sign  (positive  to  negative),  with 
z*  being  the  redshift  peak  for  an  FSRQ  with  a  luminosity  of 
1048  erg  s-1. 

The  LDDE  model  provides  a  good  fit  to  the  LAT  data  and  is 
able  to  reproduce  the  observed  distribution  in  Figure  2.  The  log- 
likelihood  ratio  test  shows  that  the  improvement  over  the  best 
PLE  model  is  significant,  with  a  chance  probability  of  ~10-6. 
Results  are  reported  in  Table  3. 

In  Figure  3,  we  subdivide  the  sample  into  four  redshift  bins 
with  comparable  number  of  sources  to  illustrate  how  the  LF 
changes.  The  evolution,  visible  as  a  shifting  of  the  peak  and  a 
change  of  the  shape  of  the  LF  between  different  bins  of  redshift, 
is  clearly  visible.  This  evolution  takes  place  mostly  below 
redshift  ~1.1  where  the  space  density  of  our  least  luminous 
objects  (i.e.,  Ly  ~  1046  erg  s-1)  increases  by  ~10  times. 
Above  this  redshift  the  variation  is  less  marked,  but  one  notices 
that 

1.  the  space  density  of  logL  =  47  objects  decreases  from 
redshift  1  to  redshift  1.5  while  that  of  log  L  =  48  FSRQs 
still  increases  (lower  left  panel,  green  versus  red  lines). 
This  indicates  that  the  space  density  of  Log  L  =  47  FSRQs 
peaks  at  a  redshift  ~  1 . 1 ;  and 


2.  likewise,  the  space  density  of  log  L  =  48  FSRQs  reaches 
its  maximum  well  below  z  =  3. 

The  best-fit  parameters  confirm  that  the  redshift  of  maximum 
space  density  increases  with  increasing  luminosity  (with  the 
power-law  index  of  the  redshift-peak  evolution  a  =  0.21  ±0.03, 
see  Equation  (16)).  This  redshift  evolution  can  be  clearly  seen  in 
Figure  4,  which  shows  the  change  in  space  density  for  different 
luminosity  classes. 

4.3.  Analysis  of  Uncertainties 

One  of  the  main  uncertainties  of  our  analysis  is  due  to  the 
incompleteness  of  the  FSRQ  sample.  In  Table  1,  there  are 
17  sources  with  associated  radio  counterparts  lacking  optical 
type  and  redshift  measurements.  A  fraction  of  these  may 
be  FSRQs.  In  addition,  there  are  29  sources  without  any 
statistically  associated  radio  counterpart.  The  lack  of  radio 
flux  means  that  these  cannot  be  FSRQ  similar  to  those  in 
the  Fermi  sample,  unless  position  errors  have  prevented  radio 
counterpart  associations.  Thus,  even  if  a  few  of  these  sources 
are  mislocalized,1 * * * * * 7  the  maximum  possible  incompleteness  of  our 
FSRQ  sample  is  on  the  order  of  20/(186+20),  i.e.,  ~10%. 

1  Fermi  counterpart  identification  depends  on  an  accurate  assessment  of 

uncertainties  in  the  source  localization.  While  the  method  properly  accounts 

for  the  catalog  uncertainties  in  the  source  position,  additional  systematic 

uncertainty  may  be  present  due  to  inaccuracies  in  the  PSF  model  or 

background  model,  particularly  in  regions  of  bright  diffuse  Galactic  emission. 

Improved  modeling  of  the  PSF  or  of  the  diffuse  emission  may  shift  the  source 

localization,  and  hence  the  counterpart  association  probabilities,  of  a  few 
sources. 
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Table  3 

Best-fit  Parameters  of  the  Luminosity  Dependent  Density  Evolution  LF 


Sample 

No.  of  Objects 

Aa 

n 

U 

Y2 

z* 

a 

Pi 

P2 

[i 

a 

ALL 

186 

3.06(±0.23)  x  104 

0.21  ±0.12 

0.84  ±  0.49 

1.58  ±0.27 

1.47  ±  0.16 

0.21  ±  0.03 

7.35  ±  1.74 

-6.51  ±  1.97 

2.44  ±  0.01 

0.18  ±0.01 

ALLb 

208 

2.82(±0.19)  x  104 

0.26  ±0.12 

0.87  ±  0.53 

1.60  ±  0.27 

1.42  ±0.15 

0.20  ±  0.03 

8.21  ±  1.78 

-5.66  ±  1.73 

2.42  ±  0.01 

0.19  ±0.01 

ALLC 

186 

8.72(±0.63)  x  103 

0.38  ±0.16 

0.89  ±  0.70 

1.60  ±0.30 

1.38  ±0.18 

0.18  ±0.03 

7.71  ±  1.84 

-4.44  ±  1.78 

Notes. 

a  In  units  of  10-13  Mpc-3  erg-1  s. 

b  22  unassociated  sources  were  included  in  this  sample  by  drawing  random  redshifts  from  the  observed  redshift  distribution  of  FSRQs. 
c  Derived  using  the  detection  efficiency  for  curved  reported  in  Figure  16. 
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Figure  3.  LF  of  the  Fermi  FSRQs  in  different  bins  of  redshift,  reconstructed  using  the  iVobs/Mndi  method.  The  lines  represent  the  best-fit  LDDE  model  of  Section  4.2. 
To  highlight  the  evolution,  the  LF  from  the  next  lower  redshift  bin  is  overplotted  (dashed  lines). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Figure  4.  Growth  and  evolution  of  different  luminosity  classes  of  FSRQs.  Note  that  the  space  density  of  the  most  luminous  FSRQs  peaks  earlier  in  the  history  of 
the  universe  while  the  bulk  of  the  population  (i.e.,  the  low  luminosity  objects)  are  more  abundant  at  later  times.  The  range  of  measured  distribution  is  determined  by 
requiring  at  least  one  source  within  the  volume  (lower  left)  and  sensitivity  limitations  of  Fermi  (upper  right). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 
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The  standard  way  to  account  for  this  incompleteness  is  to 
correct  upward  the  normalization  of  the  LF  as  to  reflect  the 
likely  real  number  of  FSRQs  (associated  or  not)  in  our  sample. 
Considering  extra  information  about  these  sources,  this  likely 
incompleteness  is  even  smaller  than  the  10%  above.  First,  we 
find  that  only  50%  of  all  the  RL  identified  sources  in  our 
sample  are  FSRQs.  This  suggests  that  only  8  of  the  17  radio 
identified  sources  are  FSRQs,  with  the  remainder  being  BL  Lac 
objects  and  lower  luminosity  AGNs.  A  similar  argument  can 
be  made  based  on  the  y- ray  spectral  index.  The  median  index 
for  FSRQs  is  T  =  2.44  while  only  9%  of  our  BL  Lac  objects 
have  such  soft  spectra.  There  are  4  such  soft  sources  in  our 
set  of  17  radio  sources.  Conservatively  assuming  that  all  are 
FSRQs,  we  infer  that  twice  this  number,  namely  eight,  FSRQs 
are  in  the  radio-detected  sample.  Of  the  29  sources  without 
radio  counterparts,  four  are  Blazar  4 ANTI- Associations”  (see 
Abdo  et  al.  2009,  2010g,  for  details),  where  we  can  definitively 
exclude  any  flat- spectrum  radio  source  bright  enough  for  a 
Fermi- type  blazar.  While  most  of  these  29  sources  lack  the  very 
deep  radio  observations  to  make  an  ANTI- Association,  all  but 
five  have  been  classified  as  pulsar  candidates,  based  on  gamma- 
ray  spectral  curvature  and  lack  of  variability.  We  thus  suspect 
that  virtually  all  of  these  sources  are  other  classes  of  gamma-ray 
emitters,  e.g.,  pulsars  yet  to  be  discovered,  starburst  galaxies, 
etc.  Our  conclusion  is  that  the  likely  incompleteness  is  only 
8/(186+8)  =  4%.  Conservatively  adding  a  few  nominally  radio¬ 
quiet  sources  from  erroneous  LAT  localizations  may  allow  5% 
incompleteness.  This  correction  has  been  applied  in  Sections  4. 1 
and  4.2. 

One  might  be  concerned  that  the  handful  of  LAT  sources  that 
could  be  unidentified  FSRQs  occupy  an  unusual  distribution  in 
redshift  space  and  could  thus  bias  the  LF.  While  this  is  possible, 
in  the  present  detected  set  there  is  no  clear  redshift  trend  with 
radio  flux.  There  is  a  dearth  of  blazars  at  high  redshift  that  are 
very  bright  in  the  optical,  but  the  optically  faintest  blazars  in  the 
set  with  spectroscopic  redshifts  happen  to  have  z  <  L  Thus,  we 
can  assume  that  the  handful  of  missing  blazars  would  show  a 
similar  redshift  distribution  to  the  identified  sources.  With  this 
assumption,  even  if  8  of  the  RL  sources  are  FSRQs  and  half  of 
the  29  sources  without  radio  counterparts  are  also  FSRQs,  we 
find  that  the  maximum  plausible  incompleteness  of  ~22  (i.e., 
8+29/2)  objects  is  insufficient  to  perturb  any  LF  parameter  out 
of  the  range  allowed  by  the  present  statistical  uncertainty  in  the 
fit  (see  Table  3).  Thus  incompleteness  is  not  a  significant  source 
of  uncertainty  in  our  study.  In  the  Appendix,  we  also  argue 
that  other  possible  sources  of  systematic  uncertainty  (detection 
efficiency,  blazar  variability,  and  absorption  by  extragalactic 
background  light,  EBL)  do  not  contribute  significantly  to  the 
uncertainty  in  the  functional  form  or  the  derived  parameters  of 
the  LF. 


4  A.  Comparison  with  Previous  Results 
4.4.1.  The  Local  Luminosity  Function 


the  maximum  allowed  volume  for  a  given  source  is  defined  as 


Vmax  — 


n(Li9z) 


e(z,  Li)  dV 

efanin,  Li)  dz 


(17) 


where  Lt  is  the  source  luminosity,  Q(L;,  z)  is  the  sky  coverage, 
Zmax  is  the  redshift  above  which  the  source  drops  out  of  the 
survey,  and  e(z,  Lt)  is  the  evolution  term  of  Equation  (15) 
normalized  (through  e(zmin ,  Lt))  at  the  redshift  zm in  to  which 
the  LF  is  to  be  de-evolved.  The  LF  de-evolved  at  Zmin  femin  =  0 
in  this  case)  is  built  using  the  standard  1  /Vmax  method  (Schmidt 
1968). 

In  order  to  gauge  the  uncertainties  that  the  different  methods 
might  introduce  in  the  determination  of  the  local  LF  we  consider 
also  an  alternate  method.  We  perform  a  Monte  Carlo  simula¬ 
tion,  drawing  1000  series  of  parameters  from  the  covariance 
matrix  of  the  best-fit  LDDE  model  described  in  Section  4.2. 
Using  the  covariance  matrix  ensures  that  parameters  are  drawn 
correctly,  taking  into  account  their  correlations.  The  re- sampled 
parameters  are  used  to  compute  the  ± la  error  of  the  LF  at  red¬ 
shift  0.  This  is  reported  in  Figure  5.  There  is  very  good  agreement 
with  the  local  LFs  using  this  method  and  the  1  /Vmax  approach. 
The  gray  band  in  Figure  5  shows  the  true  statistical  uncertainty 
on  the  space  density  that  the  1  /Vmax  method  (applied  using  the 
best-fit  parameters)  is  not  able  to  capture. 

We  find  a  local  LF  described  by  a  power  law  with  index 
of  1.6-1. 7,  for  Fioo  <  1047  erg  s-1,  steepening  at  higher 
luminosity.  Thus,  the  local  LF  can  be  parameterized  as  a  double 
power  law: 


where  A  =  (3.99  ±  0.30)  x  10"11,  L*  =  0.22  ±  0.30,  Yl  = 
1.68  zb  0.17,  Y2  =  3.15  ±  0.63,  and  both  Ly  and  L*  are  in 
units  of  1048  erg  s-1.  Other  models  (e.g.,  a  Schecter  function, 
a  simple  power  law,  etc.)  do  not  generally  provide  as  good  of 
a  fit  to  the  data.  The  values  of  the  low-luminosity  index  y\  and 
the  high-luminosity  index  y2  found  here  are  in  good  agreement 
with  those  of  1.63  ±  0.16  and  2.3  ±  0.3  reported  by  Padovani 
et  al.  (2007)  for  the  DRBXS  survey  of  FSRQs. 

Values  very  similar  to  those  found  here  were  also  reported  for 
a  radio  FSRQ  sample  by  Dunlop  &  Peacock  (1990),  who  find8 * 
Yi  =  1.83  and  y2  =  2.96.  The  Fermi  LF  low-luminosity  index 
(i.e.,  Y\)  is  flatter  than  that  determined  using  EGRET  blazars  by 
Narumoto  &  Totani  (2006)  as  is  apparent  in  Figure  5.  However,  a 
re-analysis  of  the  same  data  sets  employing  the  blazar  sequence 
(Fossati  et  al.  1998)  to  model  the  blazar  SEDs  found  a  low- 
luminosity  index  in  the  1. 8-2.1  range  (Inoue  &  Totani  2009). 
Also,  in  a  more  recent  work,  Inoue  et  al.  (2010)  modified  their 
SED  models  to  be  able  to  reproduce  TeV  data  of  known  blazars. 
Their  LF  at  redshift  0  (see  Figure  5)  is  found  to  be  in  relatively 
good  agreement  with  that  found  here  for  the  Fermi  sample. 


The  local  LF  is  the  luminosity  function  at  redshift  0.  For  an 
evolving  population,  the  local  LF  is  obtained  by  de-evolving  the 
luminosities  (or  the  densities)  according  to  the  best-fit  model. 
This  is  generally  done  using  the  1/Vmax  method  of  Schmidt 
(1968),  as  reported  for  example  by  Della  Ceca  et  al.  (2008). 
However,  since  the  best  representation  of  the  LF  is  the  LDDE 
model,  the  maximum  volume  has  to  be  weighted  by  the  density 
evolution  relative  to  the  luminosity  of  the  source.  In  this  case, 


4 A. 2.  The  Luminosity  Function  at  Redshift  1 

Figure  6  shows  the  luminosity  of  FSRQ  at  redshift  1  compared 
to  predictions  from  recent  models.  It  is  apparent  that  the 
evolution  of  the  Fermi  LF  is  stronger  than  predicted  by  any 
of  these  models.  The  increase  in  space  density  from  redshift  0 


8  Their  definition  of  local  luminosity  function  and  Equation  (18)  differ  by  a 

1/Lj,  (or  1/P  in  their  paper)  term.  Thus,  we  added  1.0  to  their  exponents. 


8 


The  Astrophysical  Journal,  751:108  (20pp),  2012  June  1 


Ajello  et  al. 


10'3  10’2  10’1  1  10 

L^l 048  erg  s1] 


Figure  5.  Local  (z  =  0)  LF  of  the  Fermi  FSRQs  as  derived  from  the  best-fit  LDDE  model  in  Section  4.2  (solid  line).  The  gray  band  represents  the  ±1  a  uncertainty 
computed  as  described  in  the  text.  The  long-  and  short-dashed  lines  show  the  LF  models  based  on  the  EGRET  blazars  derived  by  Narumoto  &  Totani  (2006)  and  Inoue 
et  al.  (2010),  respectively.  The  dash-dotted  line  shows  the  prediction  from  the  model  of  Stecker  &  Venters  (2011). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Figure  6.  LF  of  the  Fermi  FSRQs  at  redshift  1.0  as  derived  from  the  best-fit  LDDE  model  in  Section  4.2  (solid  line).  The  gray  band  represents  the  ±1  a  uncertainty 
computed  with  the  method  described  in  the  text.  The  long-dashed  and  short-dashed  lines  show  the  LF  models  based  on  the  EGRET  blazars  derived,  respectively,  by 
Narumoto  &  Totani  (2006)  and  Inoue  et  al.  (2010).  The  dash-dotted  line  shows  the  prediction  from  the  model  of  Stecker  &  Venters  (2011). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 
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to  1  for  a  source  with  a  luminosity  of  1048  erg  s-1  is  almost  a 
factor  ~150.  This  dramatic  increase  is  not  seen  in  the  evolution 
of  radio-quiet  AGNs  (e.g.,  Ueda  et  al.  2003;  Hasinger  et  al. 
2005)  whose  space  density  increases  by  a  factor  25-50  between 
redshifts  0  and  1.  The  increase  of  a  factor  ~60  seen  in  FSRQs 
detected  in  radio  (Dunlop  &  Peacock  1990)  is  still  slower  than 
that  of  Fermi  blazars.  This  explains  why  the  predictions  based 
on  LFs  derived  at  other  wavelengths  (see  Figure  6)  underpredict 
the  density  of  high-luminosity  Fermi  FSRQs  at  redshift  of  1 . 

5.  THE  SPECTRAL  ENERGY  DISTRIBUTIONS  OF  FSRQs 

We  may  use  the  0.1-300  GeV  LAT  spectra  of  our  uniform 
bright  Fermi  FSRQ  sample  along  with  the  15-200  keV  spectra 
measured  by  the  Swift  Burst  Alert  Telescope  (BAT)  to  charac¬ 
terize  the  high  energy  inverse  Compton  (IC)  sector  of  the  blazar 
SED.  Indeed,  since  Fermi-L AT  and  Swift -BAT  have  compara¬ 
ble  sensitivity  in  their  respective  bands9  and  since  the  two  bands 
cover  the  bulk  of  the  IC  component,  a  joint  study  allows  an  ac¬ 
curate  characterization  of  the  IC  spectrum  and  the  contribution 
of  FSRQs  to  the  cosmic  high-energy  background.  In  the  next 
sections,  we  describe  how  the  data  analysis  for  LAT  and  BAT 
was  performed. 


5.7.  Data  Analysis 

Our  goal  is  to  produce  SEDs  for  the  bright  Fermi  FSRQs 
that  cover  the  energy  range  15  keV-300  GeV  by  combining 
data  from  Swift /BAT  and  Fermi /LAT.  For  this  analysis,  we 
further  restrict  the  sample  of  Fermi's  FSRQs  to  F\qq  )  7  x 
10-8  photons  cm-2  s-1  (corresponding  to  an  energy  flux  of 
3.4  x  10“ 11  erg  s-1  cm-2  for  a  power  law  with  an  index  of 
2.45),  since  for  brighter  sources  Fermi  has  a  negligible  bias  in 
the  detected  spectral  indices  and  FSRQs  are  selected  uniformly 
(Abdo  et  al.  201  Of).  For  fainter  sources  hard- spectrum  objects 
(principally  BL  Lac  objects)  dominate  the  sample.  There  are  103 
bright  FSRQs  detected  by  Fermi  that  meet  these  criteria  (Abdo 
et  al.  2010g).  In  this  section,  we  describe  how  the  processing  of 
the  Fermi  data  and  that  of  Swift  data  has  been  performed. 

We  analyze  two  years  of  Fermi  data  using  version  V9r21 
of  the  science  tools.10  The  data  are  filtered,  removing  time 
periods  in  which  the  instrument  was  not  in  sky- survey  mode 
and  photons  whose  zenith  angle  is  larger  than  100°.  We  consider 
only  photons  collected  within  15°  of  the  source  position  with 
100  MeV  <  E  <  300  GeV.  We  employ  the  P7SOURCE_V6 
instrumental  response  function  and  perform  binned  likelihood 
analysis  using  the  gtlike  tool.  First,  a  likelihood  fit  using  a 
power-law  model  for  all  the  sources  in  the  region  of  interest 
is  performed  on  the  entire  energy  band  (100  MeV-300  GeV). 
The  parameters  (i.e.,  flux  and  photon  index)  of  all  the  sources 
within  3°  of  the  target  FSRQ,  along  with  the  normalization 
of  the  diffuse  model,  are  left  free.  More  distant  sources  have 
parameters  frozen  at  the  2FGL  measured  values  (Nolan  et  al. 
2012).  We  next  choose  30  logarithmically  spaced  energy  bins 
and  perform  a  binned  likelihood  in  each,  deriving  the  flux  of  the 
target  FSRQ  in  each  energy  bin.  During  this  exercise  the  flux  of 
the  FSRQ  is  allowed  to  vary  while  the  photon  index  is  fixed  at 
the  best-fit  value  found  for  the  whole  band.  All  the  neighboring 
sources  had  parameters  fixed  at  the  best-fit  values,  although 


9  A  Fermi  FSRQ  with  a  photon  flux  of  3  x  10  8  photons  cm  2  s  1  in  the 
100  MeV-100  GeV  band  and  a  power-law  spectrum  with  an  index  of  2.4  has 
an  energy  flux  of  1.5  x  10“ 11  erg  cm-2  s-1. 

10  http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ 


the  diffuse  emission  normalization  was  allowed  to  vary.  This 
analysis  provides  a  30  bin  100  MeV-300  GeV  energy  spectrum 
for  all  103  sources  in  the  bright  FSRQ  sample. 

Swift-BAI  is  a  coded-mask  telescope  that  has  conducted  a 
several-year  survey  in  the  15-200  keV  hard  X-ray  sky  (Gehrels 
et  al.  2004;  Barthelmy  et  al.  2005).  With  this  deep  exposure,  BAT 
reaches  a  sensitivity  of  ~  10-1 1  erg  cm-2  s_1  on  most  of  the  high- 
latitude  sky  (e.g.,  Tueller  et  al.  2008;  Ajello  et  al.  2008a,  2008c; 
Cusumano  et  al.  2010).  Blazars  represent  15%-20%  of  the 
extragalactic  source  population  detected  by  BAT  (Ajello  et  al. 
2009a).  We  use  years  of  BAT  data  to  extract  a  15-200  keV 
spectrum  for  all  the  FSRQs  in  the  Fermi  sample.  Spectral 
extraction  is  performed  according  to  the  recipes  presented  by 
Ajello  et  al.  (2008c)  and  discussed  in  detail  by  Ajello  et  al. 
(2009b)  and  Ajello  et  al.  (2010). 

Both  BAT  and  LAT  spectra  are  multiplied  by  4tt  Dl(z )2  (with 
Dl(z )  the  luminosity  distance  at  redshift  z)  to  transform  the  flux 
into  a  luminosity  and  shifted  by  (1  +  z)  to  transform  into  source 
rest- frame  SEDs. 

For  each  FSRQ,  we  fit  the  BAT  and  LAT  data  with  an 
empirical  model  of  the  following  form: 


•4; tDl(z)2  =  E2 


1 


e  Je,ec  m  4^ dl(z)2,  (19) 


where  y&  at  and  at  are  the  power-law  indices  in  the  BAT  and 
the  LAT  bands  and  Eb  and  Ec  are  the  break  and  the  cutoff  energy, 
respectively.  The  e~^E^Ec  term  allows  us  to  model  the  curvature 
that  is  clearly  visible  in  a  few  of  the  Fermi  spectra.  The  fit  is 
performed  only  for  E  <  20  GeV  to  avoid  possible  steeping  due 
to  the  absorption  of  y-ray  photons  by  the  EBL  (e.g.,  Stecker 
et  al.  2006;  Franceschini  et  al.  2008). 

Two  sample  spectra  are  shown  in  Figure  7.  It  is  apparent 
that  for  the  brightest  FSRQs,  BAT  and  LAT  are  efficient  in 
constraining  the  shape  of  the  IC  emission.  Even  when  the  BAT 
signal  is  not  significant,  the  upper  limit  from  BAT  still  provides 
useful  constraints  on  the  low  energy  curvature  of  the  SED.  In 
a  number  of  bright  sources  (see,  e.g.,  Figure  7)  the  highest- 
energy  data  point  in  BAT  at  >120  keV  is  seen  to  deviate  from 
the  baseline  fit.  This  deviation  is  at  present  not  significant 
(i.e.,  the  reduced  x2  of  the  baseline  fit  is  already  ~1.0), 
but  is  suggestive  of  a  second  component.  Observations  with 
INTEGRAL  extending  to  energies  >200  keV  might  ascertain 
the  nature  of  this  feature. 

Several  caveats  necessarily  apply  to  our  analysis.  First,  the 
BAT  and  LAT  observations  are  not  strictly  simultaneous.  LAT 
spectra  are  accumulated  over  two  years  while  the  BAT  data 
span  six  years.  In  principle  one  could  restrict  the  BAT  data 
to  the  period  spanned  by  the  Fermi  observations.  In  practice 
this  would  seriously  limit  the  BAT  sensitivity,  weakening 
constraints  on  most  of  the  spectra.  Second,  it  is  possible  that 
BAT  and  LAT  are  not  sampling  exactly  the  same  emission 
component.  In  particular,  BAT  might  be  dominated  by  IC 
emission  produced  by  the  synchrotron  self-Compton  (SSC; 
Maraschi  et  al.  1992)  component  while  the  LAT  may  be  more 
sensitive  to  external  Compton  (EC;  Dermer  &  Schlickeiser 
1993)  emission.  Ultimately,  detailed  SED  modeling  with  strictly 
simultaneous  data  would  be  needed  to  eliminate  these  concerns, 
and  such  work  is  well  beyond  the  scope  of  this  paper.  Bearing 
these  caveats  in  mind,  our  bright  sample  is  nearly  free  of 
selection  effects  other  than  the  hard  flux-limit  threshold  applied 
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Figure  7.  Two  BAT-LAT  spectra  of  famous  blazars  fitted  with  the  empirical 
model  described  in  Section  5. 


to  the  Fermi  data.  This  allows  a  detailed  study  of  the  average 
properties  of  the  high-energy  SED  of  FSRQs. 

5.2.  Correlation  of  Peak  Luminosity  and  the  Energy  of  the  Peak 

We  can  compare  the  IC  rest- frame  peak  luminosity  (i.e.,  the 
luminosity  at  the  peak  of  the  IC  component)  and  peak  energy 
from  the  SED  fits  to  the  FSRQs  in  our  sample.  As  shown 
in  Figure  8,  there  is  no  apparent  correlation  between  these 
quantities;  indeed  the  Kendall  test  gives  a  value  r  =  0.09 
indicating  no  significant  correlation.  Since  generally  neither 
Fermi  nor  BAT  directly  sample  the  high-energy  peak,  we  also  fit 
the  spectra  using  a  third-degree  polynomial  function  instead  of 
the  model  in  Equation  (19).  The  results  are  shown  in  the  lower 
panel  of  Figure  8  and  confirm  the  previous  findings. 

This  is  in  contrast  to  the  correlation  found  (but  see  also 
Nieppola  et  al.  2008)  between  the  luminosity  and  the  energy 
of  the  synchrotron  peak  of  blazars  (e.g.,  Ghisellini  et  al.  1998; 
Fossati  et  al.  1998).  This  might  imply  that  the  jet  parameters 
(e.g.,  Doppler  factor,  luminosity  of  the  target  photon  field,  etc.) 
do  not  depend  on  blazar  GeV  luminosity  or  redshift.  This  may  be 
understood  if  the  IC  peak  is  largely  controlled  by  EC  emission 
for  these  sources. 

5.3.  Average  SEDs 

It  is  useful  to  estimate  the  average  SED  of  FSRQs,  partic¬ 
ularly  for  estimating  the  contribution  of  FSRQs  to  the  extra- 
galactic  gamma-ray  background.  First  we  define  the  bolometric 


log(E)  [MeV] ' 


Figure  8.  Peak  luminosity  vs.  the  energy  of  the  peak  for  the  complete  sample 
of  FSRQs  discussed  in  Section  5.  The  upper  plot  shows  the  value  derived 
using  a  double  power  law  with  exponential  cutoff  while  the  lower  panel  shows 
parameters  derived  using  a  third-degree  polynomial  function. 

luminosity  as  the  luminosity  in  the  10  keV-300  GeV  band11  and 
divide  the  sources  into  four  bins  of  bolometric  luminosity  with 
approximately  the  same  number  of  objects  in  each  bin.  In  these 
luminosity  bins  we  compute  the  average  of  the  logarithm  of  the 
spectral  luminosity  at  each  energy.  Associated  errors  on  this  av¬ 
erage  spectrum  are  computed  using  the  Jackknife  technique.  In 
this  framework,  we  neglect  uncertainties  due  to  different  level 
of  the  energy  density  of  the  EBF  which  would  affect  mostly  the 
high-energy  part  of  the  SED  (i.e.,  ^20  GeV). 

Figure  9  shows  the  average  rest-frame  SED  for  the  FSRQ 
sample  in  the  four  luminosity  bins.  This  plot  confirms  the  lack 
of  a  systematic  correlation  of  the  peak  luminosity  and  energy. 
Indeed,  all  the  averaged  SEDs  show  a  peak  in  the  10-100  MeV 
band  and  their  shape  does  not  change  much  with  luminosity. 

To  transform  luminosities  between  observed  and  rest-frame 
we  need  the  k-correction,  along  with  its  redshift  variation,  shown 
in  Figure  10.  In  practice,  there  is  little  difference  between  the 
^-correction  for  the  average  SED  computed  here  (even  applying 
EBF  absorption,  e.g.,  Franceschini  et  al.  2008)  and  one  com¬ 
puted  for  a  power  law  (i.e.,  (1  +  z)T~2)  with  a  photon  index  of 

2.4.  Only  at  large  redshifts  do  the  two  k-corrections  start  to  dif¬ 
fer;  this  difference  is  only  ~5%  at  a  redshift  of  4.  We  find  that 
using  a  power-law  index  of  2.37  and  taking  into  account  EBF 
absorption  allow  us  to  reproduce  correctly  the  ^-correction  up 
to  redshift  ~6. 


11  The  best  fit  is  extrapolated  from  20  GeV  to  300  GeV. 
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Redshift-corrected  Energy  [MeV] 

Figure  9.  Average  rest-frame  spectral  energy  distributions  for  four  representa¬ 
tive  FSRQ  luminosity  classes;  see  Section  5.3.  In  each  SED,  the  band  represents 
the  1  o  uncertainty  on  the  average.  This  does  not  reflect  the  uncertainty  connected 
with  the  level  of  the  extragalactic  background  light. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


z 


Figure  10.  Ratio  of  source  rest-frame  luminosity  to  observed  luminosity  (i.e., 
k-correction)  as  a  function  of  redshift  for  the  average  SED  shape  reported  in 
Section  5.3  and  for  two  generic  power  laws. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


6.  THE  CONTRIBUTION  TO  THE  ISOTROPIC 
GAMMA-RAY  BACKGROUND 

The  nature  of  the  diffuse  gamma-ray  background  at  GeV 
energies  remains  one  of  the  most  interesting  problems  in 
astrophysics.  The  presence  of  an  isotropic  component  was  first 
indicated  by  the  OSO-3  satellite  (Kraushaar  et  al.  1972)  and 
confirmed  by  SAS-2  and  EGRET  (respectively  Fichtel  et  al. 
1975;  Sreekumar  et  al.  1998).  This  isotropic  component  is 
normally  referred  to  as  the  IGRB.  Fermi  recently  provided  a 
refined  measurement  of  this  isotropic  component  showing  that 
it  can  be  adequately  described  as  a  single  power  law  with  an 
index  of  2.4  in  the  200  MeV-100  GeV  energy  range  (Abdo 
et  al.  2010e). 

Blazars,  representing  the  most  numerous  identified  source 
population  in  the  EGRET  (Hartman  et  al.  1999)  and  Fermi 
(Ackermann  et  al.  2011)  catalogs,  are  expected  to  produce 
a  substantial  fraction  of  the  IGRB.  Typical  predictions  range 
from  20%  to  30%  (Chiang  &  Mukherjee  1998;  Miicke  &  Pohl 
2000;  Narumoto  &  Totani  2006;  Dermer  2007;  Inoue  &  Totani 
2009)  to  100%  (Stecker  &  Salamon  1996;  Stecker  &  Venters 
2011).  Analysis  of  the  source-count  distribution  showed  that  for 
Fioo  ^  10-9  photons  cm-2  s-1  the  contribution  of  unresolved 
blazars  to  the  IGRB  is  ~16%  in  the  100  MeV-100  GeV  band 
(Abdo  et  al.  201  Of).  Since  the  source  count  distributions  show  a 
strong  break  at  a  flux  of  F\  oo  ~6x  10-8  photons  cm-2  s_1,  it 
was  concluded  that  extrapolating  the  source  counts  to  zero  flux 
would  produce  ~23%  of  the  IGRB. 

Now,  with  a  measured  LF  we  can  more  robustly  evaluate  the 
emission  arising  from  faint  FSRQs.  In  addition,  the  FSRQ  SED 
shape  study  of  the  previous  section  also  allows  improvement 
over  the  simple  power-law  type  spectra  assumed  by  Abdo  et  al. 
(201  Of). 

The  contribution  of  “unresolved”  FSRQs  to  the  IGRB  can  be 
estimated  as 


Tegb 
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where  the  limits  of  integration  are  the  same  as  those  of 
Equation  (5)  and  Fy(LY ,  z)  is  the  flux  of  a  source  with  rest-frame 
luminosity  Ly  at  redshift  z.  Since  we  are  interested  in  the  diffuse 
flux  not  yet  resolved  by  Fermi  (Abdo  et  al.  2010e)  the  term 
(1— Q(r,  EK)/Qmax)  takes  into  account  the  photon  index  and 
source  flux  dependence  of  the  LAT  source  detection  threshold 
(see  Abdo  et  al.  201  Of  for  more  details). 

In  Equation  (20)  setting  Q(r,  Fy)/£lmax  =  0  allows  us  to 
compute  the  total  y- ray  emission  arising  from  the  FSRQ  class. 
The  result  is  3.13^q3275  x  10-6  photons  cm-2  s-1  sr-1  in  the 
100  MeV-100  GeV  band.  This  value  should  be  compared  with 
the  total  sky  intensity12  of  ~1.4  x  10-5  photons  cm-2  s-1  sr-1, 
which  includes  IGRB  plus  detected  sources13  (Abdo  et  al. 
2010e).  Thus,  FSRQs  make  21.7+2,57%  of  this  total  intensity. 

If  one  considers  the  contribution  only  from  the  FSRQs  that 
Fermi  has  not  detected14  (i.e.,  unresolved  sources)  then  this  be¬ 
comes  9.66+\ qq  x  10-7  photons  cm-2  s_1  sr-1  (with  a  maximum 
systematic  uncertainty  of  3  x  10-7  photons  cm-2  s_1  sr-1;  see 
the  Appendix).  This  represents  9.3^\60%  of  the  IGRB  intensity 
in  the  0.1-100  GeV  band  (Abdo  et  al.  2010e).  From  above  it 
is  also  clear  that  Fermi  has  already  resolved  more  than  50% 
of  the  total  flux  arising  from  the  FSRQ  class.  Figure  1 1  shows 
this  contribution.  The  possible  presence  of  EC  components  in 
the  SEDs  of  FSRQs  makes  the  estimate  between  200  keV  and 
100  MeV  uncertain  (see  Section  5.3).  Future  observations  with 
both  Fermi  above  20  MeV  and  INTEGRAL  above  200  keV  and 
physical  modeling  of  blazar  spectra  might  substantially  reduce 
this  uncertainty. 

Even  the  (disfavored)  PEE  model  cannot  accommodate  a 
much  larger  contribution  of  FSRQs  to  the  IGRB.  Indeed,  in 
this  case  the  contribution  of  unresolved  FSRQs  would  be 
1.37  x  10-6  photons  cm-2  s_1  sr-1  (or  ~13%  of  the  IGRB 
intensity).  In  any  case,  the  sharp  falloff  of  the  FSRQ  contribution 
at  high  and  low  energies  (Figure  11)  requires  that  other  sources, 


12  This  does  not  include  the  Galactic  diffuse  emission. 

13  This  flux  is  obtained  summing  the  IGRB  and  the  resolved-source  fluxes 
provided  in  Table  1  of  Abdo  et  al.  (2010e)  and  extrapolating  their  sum  to 

100  MeV. 

14  This  is  important  for  the  comparison  of  the  unresolved  emission  of  FSRQs 
with  the  measurement  of  the  IGRB  which  does  not  include  point  sources 
detected  by  Fermi  during  the  first  ~9  months  of  observations  (Abdo  et  al. 
2010e). 
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Energy  [MeV] 

Figure  11.  Contribution  of  unresolved  (top)  and  total  (resolved  plus  unresolved,  bottom)  FSRQs  to  the  diffuse  extragalactic  background  (blue  line)  as  determined 
by  integrating  the  luminosity  function  coupled  to  the  SED  model  derived  in  Section  5.3.  The  hatched  band  around  the  best-fit  prediction  shows  the  lcr  statistical 
uncertainty  while  the  gray  band  represents  the  systematic  uncertainty. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


e.g.,  BL  Lac  objects  and  starburst  galaxies  make  significant 
contributions  to  the  IGRB  intensity. 

7.  BEAMING:  THE  INTRINSIC  LUMINOSITY  FUNCTION 
AND  THE  PARENT  POPULATION 

The  luminosities  L  defined  in  this  work  are  apparent  isotropic 
luminosities.  Since  the  jet  material  is  moving  at  relativistic  speed 
(y  >1),  the  observed,  Doppler  boosted,  luminosities  are  related 
to  the  intrinsic  values  by 

L  =  8P  X  (21) 


where  Jjf  is  the  intrinsic  (unbeamed)  luminosity  and  8  is  the 
kinematic  Doppler  factor 

8  =  (y-  A2  -  lcos  0)~\  (22) 

where  y  —  (1  —  /32)-1/2  is  the  Lorentz  factor  and  fi  =  v/c  is  the 
velocity  of  the  emitting  plasma.  Assuming  that  the  sources  have 
a  Lorentz  factor  y  in  the  y\  ^  y  ^  72  range  then  the  minimum 
Doppler  factor  is  8m jn  =  y2_1  (when  0  =  90°)  and  the  maximum 
is  <Smax  =  72  t \J — 1  (when  0  =0°).  We  adopt  a  value  of  p  =  4 
that  applies  to  the  case  of  jet  emission  from  a  relativistic  blob 
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radiating  isotropically  in  the  fluid  frame.  The  case  p  =  4  applies 
also  to  spherical  blobs  if  the  observed  emission  is  dominated  by 
the  SSC  component,  while  a  value  of  p  =  5-6  should  be  adopted 
if  the  emission  is  due  to  EC  (Dermer  1995).  However,  as  shown 
later  these  values  imply  extremely  small  isotropic  rest-frame 
luminosities. 

Beaming  is  known  to  alter  the  shape  of  the  intrinsic  LF.  Urry 
&  Shafer  (1984)  provide  an  analytic  solution  to  the  case  where 
the  intrinsic  LF  is  a  single  power  law  and  the  jets  have  a  single 
Lorentz  factor.  In  Urry  &  Padovani  (1991)  the  intrinsic  LF  may 
be  a  double  power  law  and  a  distribution  of  Lorentz  factor  is 
considered. 

In  this  section,  we  will  determine  the  intrinsic  LF  of  the  Fermi 
FSRQs  and  their  Lorentz  and  Doppler  factor  distributions.  In 
what  follows  we  adopt  the  formalism  and  symbols  of  Lister 
(2003)  and  Cara  &  Lister  (2008). 

We  begin  by  defining  the  intrinsic  LF  as 

O  OS?)  =  kiJzT5  (23) 

valid  in  the  ^  ^  range.  The  probability  of  observing 

a  beamed  luminosity  L  given  a  Doppler  factor  8  is  (see  also 
Lister  2003) 

P(L,  8)  =  P§(8)  *  0(«S?)— — ,  (24) 

dL 

where  Ps(8)  is  the  probability  density  for  the  Doppler  8. 
Assuming  a  random  distribution  for  the  jet  angles  (i.e.,  Pq  = 
sin  0 ),  this  results  in 

Ps(S)  =  f  Pr(Y)Pe(0) 

From  here  it  follows  that 

9  fY2  Py(y) 

Ps(S)  =  S~2  /  dy,  (26) 

Jm  v  y1  -  i 

where  PY(y)  is  the  probability  density  for  y  and  the  lower  limit 
of  integration  f(8)  depends  on  the  Doppler  factor  value  and  is 
reported  in  Equation  (A6)  in  Lister  (2003).  Integrating  over  8 
yields  the  observed  LF  of  the  Doppler  beamed  FSRQs: 

fS2(L) 

<&(L)  =  klL~B  \  Ps(S)SpiB-l\  (27) 

Jsi(L) 

where,  as  in  Cara  &  Lister  (2008),  the  limits  of  integration  are 
8i(L)  =  min{<5max,  max(Smin,  (L/Sf2)1,p)}  (28) 

82(L)  =  max{(5min,  min(<W,  (L/Jzfi)1^)}-  (29) 

In  this  way,  by  fitting  Equation  (27)  to  the  Fermi  Doppler  boosted 
LF,  it  is  possible  to  determine  the  parameters  of  the  intrinsic  LF 
and  of  the  Lorentz-factor  distribution. 

We  assume  that  the  probability  density  distribution  for  y  is  a 
power  law  of  the  form 

Py(Y)  =  Cy\  (30) 

where  C  is  a  normalization  constant  and  the  function  is  valid 
for  yi  ^  y  ^  ]/2-  Later  we  will  discuss  other  parameterizations 
of  the  distribution  of  Lorentz  factors.  Here,  we  assume  y\  =  5 
and  72  =  40  as  this  is  the  range  of  Lorentz  factors  observed 


Table  4 

Parameters  of  the  Beaming  Models  Described  in  the  Text 


Parameter 

Value 

Value 

k 

-2.03  ±  0.70 

-2.43  ±0.11 

k\ 

5.1  ±  0.5a 

5.0  ±  0.5b 

B 

3.04  ±  0.08 

3.00  ±  0.08 

Yl 

5 

5 

Y2 

40 

40 

& i 

1040 

1038 

1044 

1042 

P 

4 

5 

X2/dof 

1.3 

1.5 

Notes.  Parameters  without  an  error  estimate  were  kept  fixed 
during  the  fitting  stage. 
a  In  units  of  10-23. 
b  In  units  of  10-26. 

for  y  -loud  blazars  (e.g.,  Lahteenmaki  &  Valtaoja  2003;  Lister 
et  al.  2009b;  Savolainen  et  al.  2010).  While  the  largest  intrinsic 
luminosity  (£2)  can  be  set  free,  the  lowest  one  depends  on  the 
value  of  p  chosen,  i.e.,  from  Equation  (21)  the  beaming  factor 
defines  the  intrinsic  luminosity  corresponding  to  the  apparent 
isotropic  luminosity  we  observe.  For  p  —  4  and  p  —  5,  C\  has 
to  be  set  equal  to  1040  erg  s-1  and  1038  erg  s-1,  respectively.  We 
set  £2  =  1 04  £  1 ,  but  this  choice  has  hardly  any  impact  on  the 
results. 

The  free  parameters  of  the  problem  are  the  normalization  (k\) 
and  the  slope  ( B )  of  the  intrinsic  LF  and  the  slope  k  of  the  Lorentz 
factor  distribution.  We  have  fitted  Equation  (27)  to  the  Fermi  LF 
de-evolved  at  redshift  0  derived  in  Section  4.4.1.  Figure  12 
shows  how  the  best-fit  beaming  model  reproduces  the  local  LF 
of  FSRQs  measured  by  Fermi.  From  the  best  fit  we  derive,  for 
the  p  =  4  case,  an  intrinsic  LF  slope  of  B  =  3.04  zb  0.08  and  an 
index  of  the  Lorentz-factor  distribution  of  k\  =  —2.03  ±  0.70. 
The  fit  values  are  summarized  in  Table  4.  The  Lorentz-factor 
distribution  implies  an  average  Lorentz  factor  (defined  by  the 
range  of  Lorentz  factors  given  in  Table  4)  for  Fermi  blazars  of 
y  =  11.7+22,  reasonable  agreement  with  measured  values 
(see,  e.g.,  Ghisellini  et  al.  2009).  The  index  of  the  Lorentz-factor 
distribution  is  in  agreement  with  k\  — — 1.5  found  for  the  CJ-F 
survey  (Lister  &  Marscher  1997).  The  parameters  for  the  p  =  5 
case  are  very  similar  to  those  of  the  p  =  4  case,  but  the  reduced 
X2  is  slightly  worse  (see  Table  4).  Nevertheless,  the  predictions 
of  the  two  models  are  in  agreement  and  we  find  again  that 
the  average  bulk  Lorentz  factor15  is  y  —  10.2+2  4*  As  noted 
already  the  extreme  Doppler  boosting  (85)  requires  the  intrinsic 
luminosities  to  be  small,  i.e.,  1038  erg  s-1  <  ££  ^  1042  erg  s-1. 

From  the  ratio  of  the  integrals  of  the  beamed  and  intrinsic 
LF  we  derive  that  the  Fermi  FSRQs  represent  only  ~0.1%  of 
the  parent  population.  The  average  space  density  of  LAT  FSRQs 
(derived  from  the  LF,  Section  4.2)  is  1 .4  Gpc-3 ,  implying  that  the 
average  space  density  of  the  parent  population  is  —1500  Gpc-3. 
Our  model  also  allows  us  to  determine  the  distribution  of  jet 
angles  with  respect  to  our  line  of  sight.  This  is  found  to  peak 
at  — 1?0  (Figure  13).  While  FSRQs  can  still  be  detected  at 
large  (i.e.,  >10°)  off-axis  angles  for  reasonably  low  y  factors 
(—5-7),  most  FSRQs  detected  by  Fermi  are  seen  at  angles 
less  than  5°  from  the  jet  axis.  Owing  to  their  larger  space 
density  (see  Figure  12),  misaligned  jets  produce  a  non-negligible 
diffuse  emission.  From  our  model  it  is  found  that  the  ratio 


15  We  remind  the  reader  that  the  distribution  of  Lorentz  factors  is  limited  to  be 
in  the  5-40  range  (see,  i.e.,  Table  4). 


dy 


/ 


as* 
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Figure  12.  Fermi's  LF  de-evolved  at  redshift  0  and  the  best-fit  beaming  model  (for  p  =  4;  see  the  text)  described  in  Section  7.  The  red  continuous  line  shows  the 
predicted  space  density  of  misaligned  jets. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


between  the  diffuse  emission  contribution  of  misaligned  jets 
and  that  of  blazars  (at  redshift  0)  is  ~30%.  This  has  obvious 
consequences  for  the  generation  of  the  IGRB.  In  fact  nearly 
all  of  the  flux  produced  by  radio  galaxies  is  unresolved,  with 
only  two  steep- spectrum  radio  quasars,  and  two  FR  II  and 
seven  FR  I  radio  galaxies  detected  with  Fermi-LAI  (Abdo 
et  al.  2010b).  All  the  results  reported  above  apply  to  both  the 
p  =  4  and  p  =  5  models.  Finally  we  also  tested  different 
parameterizations  of  the  distribution  of  Lorentz  factors,  i.e., 
Equation  (30).  We  used  a  linear,  an  exponential,  and  a  Gaussian 
distribution.  The  only  acceptable  fit  (but  statistically  worse  than 
the  power  law)  is  produced  by  an  exponential  distribution  of 
the  form:  P(y)  oc  e~omy  which  is  similar  to  our  power-law 
model  for  y  >  10.  We  thus  conclude  that  the  parameterization 
of  the  Lorentz  factor  distribution  with  a  power-law  model  is  a 
reasonable  assumption. 

8.  DISCUSSION  AND  CONCLUSIONS 

In  this  paper,  we  examine  the  properties  of  y  -ray-selected 
FSRQs  using  data  from  the  Fermi-LAI  instrument.  Our  work 
relies  on  a  nearly  complete,  flux-limited  sample  of  1 86  FSRQs 
detected  by  Fermi  at  high  significance  and  large  Galactic  latitude 
during  the  first  year  of  operation.  This  analysis  explores  several 
of  the  properties  of  FSRQs;  here  we  discuss  and  summarize  our 
findings. 

8.1.  Beamed  Luminosity  Function 

The  redshift-zero  LF  of  Fermi  FSRQs  is  well  described  by 
a  double  power-law  model,  typical  for  the  LF  of  both  radio¬ 
quiet  and  RL  AGNs.  At  mid-to-high  luminosities  there  is  good 
agreement  between  the  Fermi  LF  and  that  determined  using 
EGRET  data  (e.g.,  Narumoto  &  Totani  2006;  Inoue  et  al.  2010). 
At  luminosities  <1046  erg  s-1  the  FSRQ  LF  appears  to  be 
slightly  flatter  than  in  previous  studies. 


Figure  13.  Distribution  of  viewing  angles  with  respect  to  the  jet  axis  for  Fermi 
FSRQs. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


The  space  density  of  LAT-detected  FSRQs  increases  dramat¬ 
ically  with  redshift,  growing  by  50-100  times  by  z  =  L5. 
Describing  the  evolution  of  the  LF  as  simple  luminosity  evolu¬ 
tion  (PLE  model),  there  are  strong  indications  that  the  evolution 
must  cut  off  for  z  ^  1.6.  After  this  redshift,  the  space  density 
of  blazars  starts  to  decrease  quickly. 

A  simple  PLE  does  not  fully  explain  the  Fermi  data.  In  par¬ 
ticular  the  source  count  distribution  is  not  well  modeled.  Since 
there  is  evidence  that  low-  and  high-luminosity  sources  have  dif¬ 
ferent  redshift  peaks,  we  consider  a  more  sophisticated  model 
where  the  evolution  is  primarily  in  density,  but  objects  with  dif¬ 
ferent  luminosity  are  allowed  to  have  different  redshift  peaks. 
This  so-called  LDDE  model  explains  well  the  evolutionary  be¬ 
havior  of  (radio-quiet)  AGNs  selected  in  the  X-ray  band  (Ueda 
et  al.  2003;  Hasinger  et  al.  2005)  and  was  also  suggested  by 
Narumoto  &  Totani  (2006)  to  describe  the  LF  of  EGRET  blazars. 
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Redshift 


Figure  14.  Luminosity  density  as  a  function  of  redshift  produced  by  the  Fermi 
FSRQs.  The  gray  band  represents  the  1  o  statistical  uncertainty  around  the  best- 
fit  LF  model. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Figure  15.  Number  density  of  LAT-detected  FSRQs  as  a  function  of  redshift. 
The  gray  band  represents  the  1  a  statistical  uncertainty  around  the  best-fit  LF 
model. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


The  LDDE  model  provides  a  good  description  of  the  LF  of  the 
Fermi  FSRQs.  We  find  that  the  predictions  reported  in  the  liter¬ 
ature  (e.g.,  Narumoto  &  Totani  2006;  Inoue  et  al.  2010;  Stecker 
&  Venters  2011)  are  not  in  agreement  with  the  LF  of  Fermi 
FSRQs  at  redshift  unity.  This  is  due  to  the  fact  that  the  Fermi 
FSRQ  LF  is  found  to  evolve  more  quickly  than  the  LFs  of 
X-ray- selected  AGNs  or  radio- selected  FSRQs.  Indeed,  the 
space  density  of  Fermi  FSRQs  increases  by  a  factor  ~150  be¬ 
tween  redshifts  0  and  1  while  the  density  increase  is  at  most  a 
factor  ~60  for  the  models  discussed  above. 

The  LDDE  model  implies  that  sources  with  a  luminosity 
of  1046  erg  s_1,  1047  erg  s_1,  and  1048  erg  s-1  reach  their 
maximum  space  density  at  a  redshift  of  ~0.6,  ~0.9,  and  ~1.5, 
respectively.  It  is  clear,  then,  that  the  most  luminous  objects, 
while  lower  in  numbers,  appear  before  the  bulk  of  the  (low- 
luminosity)  population.  This  downsizing  in  luminosity,  where 
the  most  luminous  objects  are  found  at  earlier  times,  is  common 
to  all  classes  of  AGNs  (see,  e.g.,  Cowie  et  al.  1999;  Hasinger 
et  al.  2005,  and  references  therein),  but  is  observed  here  for 
the  first  time  at  gamma-rays.  This  type  of  downsizing  does  not 
necessarily  reflect  more  conventional  downsizing  in  terms  of 
increased  star  formation  activity  in  less  massive  galaxies  at  late 
times  (Cowie  et  al.  1996)  because  the  host  galaxy  and  black 
hole  masses  are  not  known,  but  corresponds  to  the  downsizing 
behavior  observed  in  submillimeter  (sub-mm)  galaxies  (Wall 
et  al.  2008)  even  though  the  underlying  astrophysics  is  different 
(on  the  one  hand,  concerning  relativistic  jets,  and,  on  the 
other,  dust-enshrouded  sub-mm  galaxies).  The  disappearance 
of  quasar-like  objects  at  late  times  might  indicate  that  accretion 
efficiency  evolves  as  a  function  of  cosmic  time  (e.g.,  Merloni 
2004).  Sanders  et  al.  (1988)  and  Di  Matteo  et  al.  (2005)  propose 
that  the  merging  of  two  massive  galaxies  leads  to,  in  addition 
to  strong  star  formation  activity,  a  burst  of  inflow  feeding  gas  to 
the  SMBH  and  initiates  a  “quasar-like”  phase.  Eventually  the 
energy  released  by  the  AGN  in  the  form  of  powerful  winds 
expels  the  gas,  quenches  star  formation,  and  starves  the  AGN. 
This  picture,  coupled  with  the  fact  that  major  mergers  become 
increasingly  rare  at  low  redshift  (e.g.,  Fakhouri  et  al.  2010; 
Kulkarni  &  Loeb  2012),  may  explain  why  quasars  are  rare  in 
the  local  universe. 

Figure  14  shows  the  energy  density  injected  in  the  universe 
(e.g.,  the  luminosity  density)  by  FSRQs  as  function  of  redshift. 
This  shows  a  broad  peak  between  a  redshift  of  1  and  2.  A  similar 
behavior  is  shown  by  the  cosmic  star  formation  history  (e.g., 


Hopkins  &  Beacom  2006)  which  peaks  around  redshifts  1-2. 
This  represents  a  strong  link  between  the  host  and  the  nucleus. 
A  noteworthy  fact  is  the  mild  evidence  for  a  fast  decline  in  the 
space  density  of  FSRQs  after  the  redshift  peak  (see  parameter 
P2  in  Table  3  and  also  Figures  4  and  15).  The  decline  seems  to 
be  as  dramatic  as  the  increase  in  space  density  leading  up  to  the 
redshift  peak.  For  comparison,  X-ray- selected  samples  of  AGNs 
show  a  much  milder  decline  ( p2  ~  —1.5)  after  the  redshift  peak 
(e.g.,  Ueda  et  al.  2003;  Hasinger  et  al.  2005;  Aird  et  al.  2010). 
However,  Schmidt  et  al.  (1995)  and  more  recently  Silverman 
et  al.  (2008)  reported  evidence  for  a  similarly  dramatic  decrease 
in  the  space  density  of  AGNs. 

One  factor  contributing  to  this  phenomenon  is  the  difficulty 
for  Fermi  to  detect  soft  sources  (Abdo  et  al.  2010f).  At  redshift 
^3  the  SED  peak  should  move  well  below  the  current  LAT 
energy  band,  making  it  difficult  to  probe  a  population  of 
extremely  soft  sources.  Increasing  the  effective  area  at  or  below 
100  MeV  may  help  uncover  such  a  population.  Because  the 
rising  part  of  the  IC  peak  is  in  the  hard  (>10  keV)  X-ray  band, 
high-redshift  objects  are  more  easily  selected  in  this  band  (see, 
e.g.,  the  Swift/BAI  results  in  Ajello  et  al.  2009a).  In  this  case 
another  strategy  would  be  to  build  a  bolometric  LF  that  uses 
both  the  y-ray  and  the  X-ray-selected  samples. 

8.2.  The  Intrinsic  Luminosity  Function 

Doppler  boosting  allows  Fermi  to  detect  many  blazars  when 
their  jet  emission  is  within  a  few  degrees  from  the  line  of  sight. 
As  shown  first  by  Urry  &  Shafer  (1984),  Doppler  boosting  is 
known  to  alter  the  shape  of  the  LF.  In  this  paper,  we  adopted 
a  formalism  that  allowed  us  to  recover  the  intrinsic  de-beamed 
LF  and  to  determine  the  distribution  of  bulk  Lorentz  factors  for 
the  Fermi  FSRQs. 

We  found  that  the  intrinsic  LF  can  be  modeled  by  a  single 
steep  power  law  with  an  index  of  3.04  zb  0.08  in  the  range  of 
intrinsic  luminosities  1040  erg  s-1  <  <  1044  erg  s_1.  The 

break  seen  in  the  beamed  LF  at  redshift  0  is  thus  produced  by 
Doppler  boosting.  The  data  cannot  be  explained  by  a  single, 
averaged,  Lorentz  factor,  but  require  a  distribution  of  Lorentz 
factors.  This  distribution  is  found  to  be  compatible  with  a  power 
law  with  an  index  of  —2.03  ±  0.70  in  the  5  ^  y  ^  40  range.  This 
yields  the  result  that  the  average  FSRQ  bulk  Lorentz  factor  is 
y  =  11 .7+22,  g°°d  agreement  with  several  studies  (Ghisellini 
et  al.  2009).  Our  model  is  able  to  predict  the  distribution  of 


16 


The  Astrophysical  Journal,  751:108  (20pp),  2012  June  1 


Ajello  et  al. 


viewing  angles  with  respect  to  the  jet  axes  of  Fermi  FSRQs. 
It  is  found,  see  Figure  13,  that  on  average  FSRQs  are  seen 
within  an  average  angle  of  ~2?3  from  the  jet  and  that  most 
are  seen  within  5°-6°.  A  few  Fermi  FSRQ  detections  may  be 
up  to  15°  off-axis  (if  these  have  low  Doppler  factors).  Fermi- 
detected  FSRQs  represent  only  ~0.1%  of  the  parent  population 
for  randomly  pointed  jets. 

Monitoring  observations  with  the  Very  Long  Baseline  Array 
(VLB A)  established  that  LAT-detected  blazars  have,  on  average, 
significantly  faster  apparent  jet  speeds  than  non-LAT  detected 
blazars  (Lister  et  al.  2009a;  Savolainen  et  al.  2010).  Their 
distribution  of  Lorentz  factors  is  in  good  agreement  with  the 
results  of  our  analysis.  Moreover,  they  report  the  distribution 
of  viewing  angles  with  respect  to  the  jet  axis  for  FSRQs 
detected  by  LAT.  From  their  study  the  average  viewing  angle  is 
2?9  =b  0?3  and  all  the  FSRQs  in  their  sample  have  0  <  5°.  There 
is  thus  substantial  agreement  between  the  VLBA  monitoring 
observations  and  the  results  of  our  beaming  model  applied  to 
y  -ray-only  data. 

The  space  density  of  FR-II  radio  galaxies  (i.e.,  the  putative 
parent  population  of  FSRQs)  is  reported  to  be  ~1580  Gpc-3 
(at  15  GHz)  and  ~2200  Gpc-3  (at  1.4  GHz)  by  Cara  &  Lister 
(2008)  and  Gendre  et  al.  (2010),  respectively.  From  our  study 
we  derive  a  space  density  of  FSRQ  parents  of  ~1500  Gpc-3  in 
substantial  agreement  with  the  numbers  above. 

Future  work  may  test  whether  the  intrinsic  properties  of 
blazars  (i.e.,  Lorentz  factor,  luminosity,  etc.)  evolve  with 
redshift.  This  will  likely  require  larger  complete  samples  and 
improved  modeling  of  selection  effects. 

8.3.  Spectral  Energy  Distribution 

Blazar  SEDs  are  characterized  by  the  typical  “two  hump” 
spectrum  where  the  low-energy  peak  is  produced  by  electrons 
radiating  via  synchrotron  and  the  high  energy  peak  is  produced 
via  IC  scattering  off  the  same  synchrotron  photons  (SSC 
scenario;  Maraschi  et  al.  1992)  and/or  external  seed  photons 
(EC  scenario;  Dermer  &  Schlickeiser  1993). 

In  this  work,  we  have  combined  quasi- simultaneous 
Swift /BAT  and  Fermi /LAT  data  to  investigate  the  empirical 
properties  of  the  IC  component  of  the  SEDs  of  the  FSRQs 
detected  by  Fermi.  All  the  SEDs  show  apparent  curvature  and 
have  a  peak  somewhere  in  the  10  MeV-1  GeV  band.  For  the 
FSRQs  in  our  sample,  we  have  determined  the  peak  luminosity 
of  the  IC  gamma-ray  component  and  the  rest-frame  peak  fre¬ 
quency  (or  peak  energy)  at  which  the  IC  component  reaches  its 
peak  luminosity.  No  correlation  is  found  either  between  peak 
luminosities  and  peak  energies  or  between  bolometric  lumi¬ 
nosities  of  the  IC  component  and  peak  energies,  as  shown  in 
Figure  8.  The  existence  of  such  correlation  has  been  claimed  in 
the  past  for  the  luminosity  and  the  energy  of  the  synchrotron 
peak  (Ghisellini  et  al.  1998;  Fossati  et  al.  1998)  for  a  sample 
of  blazars  (i.e.,  FSRQs  and  BL  Lac  objects).  Thus,  it  might  be 
that  this  correlation  (if  real)  exists  only  when  the  two  families 
of  blazars  are  joined  together  and  that  any  correlation  for  the 
FSRQs  class  is  washed  away  by  the  presence  of  the  additional 
EC  component.  Also  the  lack  of  correlation  of  the  IC  peak  lu¬ 
minosity  and  frequency  reveals  that  FSRQs  are,  unlike  BL  Lac 
objects,  part  of  a  population  with  homogenous  properties. 

We  built  average  redshift-corrected  SEDs  in  four  different 
luminosity  bins.  The  average  SEDs  are  surprisingly  similar  as 
a  function  of  luminosity  (and  redshift)  as  Figure  9  testifies. 
Approximating  the  SED  with  a  power  law  with  an  index  2.4, 
while  not  producing  the  correct  shape,  allows  the  reader  to 


compute  a  k-correction  useful  up  to  redshift  ~2.  Beyond  that 
this  approximation  is  not  valid. 

8.4.  The  Contribution  to  the  Diffuse  Background 

This  work  has  important  consequences  for  our  understand¬ 
ing  of  the  origins  of  the  diffuse  background.  As  pointed  out 
by  several  authors  (e.g.,  Inoue  &  Totani  2009)  and  determined 
in  this  work,  the  spectrum  of  the  diffuse  emission  arising  from 
FSRQs  shows  curvature,  due  to  the  curved  SEDs  of  these  ob¬ 
jects.  We  couple  our  model  SED  to  our  LF  to  predict  the  FSRQ 
contribution  to  the  10  keV  to  100  GeV  diffuse  background. 
FSRQs  produce  ~10%  of  the  cosmic  diffuse  emission  in  the 
1  MeV-1 0  GeV  band. 

Because  of  its  good  sensitivity  Fermi  has  already  re¬ 
solved  as  much  as  50%  of  the  total  flux  from  FSRQs  in  the 
100  MeV-1 00  GeV  band.  Our  analysis  indicates  that  the  contri¬ 
bution  of  unresolved  FSRQs  to  the  IGRB  (Abdo  et  al.  2010e)  is 
9.66t11'6079  x  10-7  photons  cm-2  s-1  sr-1  and  thus  only  9.3f\60% 
(±3%  systematic  uncertainty)  of  the  intensity  of  the  IGRB.  This 
analysis  is  in  good  agreement  with  the  results  reported  by  Abdo 
et  al.  (201  Of)  except  above  10  GeV  where  the  use  of  a  simple 
power  law  for  the  spectra  of  FSRQs  was  inadequate. 

Our  results  appear  in  reasonably  good  agreement  with  those 
of  Inoue  et  al.  (2010)  and  of  Inoue  &  Totani  (201 1),  particularly 
in  terms  of  spectral  shape  of  the  diffuse  emission  arising  from 
FSRQs.  These  authors  rely  on  the  EGRET  sample,  which 
contains  both  FSRQs  and  BL  Lac  objects,  so  it  is  not  surprising 
that  their  estimates  of  the  blazar  contribution  to  the  EGRB  are 
larger,  by  a  factor  ~2  at  1  GeV,  than  ours.  Finally,  our  estimate 
reported  above  is  in  good  agreement  with  the  results  of  Dermer 
(2007)  that  predicted  that  FSRQs  would  produce  ~10%-15% 
of  the  y- ray  background. 

LAT-detected  FSRQs  represent  only  ~0.1%  of  the  parent 
population  (see  Section  7)  and  thus  it  is  reasonable  to  expect  that 
misaligned  jets,  although  less  luminous  but  more  numerous,  give 
a  non-negligible  contribution  to  the  diffuse  background.  Our 
beaming  model  allowed  us  to  explore  this  scenario.  It  is  found 
that  misaligned  relativistic  jets  contribute  ~30%  of  the  diffuse 
flux  from  the  FSRQs  class  at  redshift  0.  If  the  Lorentz  factor 
distribution  does  not  change  with  redshift  then  the  contribution 
of  unresolved  FSRQs  and  their  misaligned  siblings  might  be 
around  ~2.0  x  10-6  photons  cm-2  s-1  and  thus  ~20%  of  the 
IGRB.  Recently,  Inoue  &  Totani  (2011)  predicted  that  radio 
galaxies  of  both  the  FR-I  and  FR-II  types  might  be  able  to 
account  for  ~25%  of  the  intensity  of  the  IGRB.  In  our  work 
we  found  that  FR-II  alone  could  in  principle  (see  above  caveat) 
produce  ~10%  of  the  IGRB.  It  can  be  envisaged  that  once  also 
the  contribution  to  the  IGRB  of  BL  Lac  objects  and  their  parents 
will  be  established,  the  total  y- ray  emission  from  relativistic 
jets  might  account  for  some  ~40%-50%  of  the  intensity  of  the 
IGRB.  Moreover,  star- forming  galaxies  at  lower  energies  are 
likely  to  remedy  some  of  the  differences  between  the  intensity  of 
the  IGRB  and  the  contribution  from  FSRQs  shown  in  Figure  11, 
though  contributions  from  other  source  classes  may  be  required 
to  explain  the  entire  IGRB  spectrum. 
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APPENDIX 

SYSTEMATIC  UNCERTAINTIES 

The  sources  of  systematic  uncertainties  in  this  analysis  are 
incompleteness  (i.e.,  missing  redshifts),  detection  efficiency, 
blazar  variability,  and  EBL.  A  detailed  discussion  of  some 
of  these  problems  was  already  given  in  Abdo  et  al.  (201  Of). 
Incompleteness  in  our  sample  is  very  small  and  introduces  no 
appreciable  systematic  uncertainty  as  shown  in  Section  4.3. 

A.l.  Detection  Efficiency 

The  detection  efficiency  used  to  determine  the  sky  area 
surveyed  by  Fermi  at  any  given  flux  is  very  important  in  this 
analysis  (see  Abdo  et  al.  201  Of,  for  a  detailed  discussion).  The 
detection  efficiency  used  in  this  work  was  derived  in  Abdo  et  al. 
(201  Of)  under  the  assumption  that  the  blazar  spectra  can  be 
approximated  by  a  power  law.  While  this  might  be  true  over 
a  small  energy  band,  it  becomes  a  problematic  assumption 
over  the  full  100  MeV-100  GeV  band  covered  by  LAT.  In 
this  Appendix,  we  estimate  directly  the  systematic  uncertainties 
connected  to  this  hypothesis.  We  performed  three  end-to-end 
Monte  Carlo  simulations  of  the  Fermi  sky  (see  Abdo  et  al. 
201  Of  for  details),  assigning  randomly  a  curved  spectrum  to 
each  source.  These  spectra  are  extracted  from  a  library  created 
using  the  ~100  observed  spectra  derived  in  Section  5  varying 
the  parameters  of  the  measured  spectra  within  their  errors.  The 
simulations  were  then  analyzed  to  derive  the  detection  efficiency 
reported  in  Figure  16.  In  particular  in  order  to  detect  a  source,  an 
ML  fit  with  a  power-law  spectrum  is  performed.  This  is  done  in 
order  to  reproduce  the  inherent  systematic  uncertainty  of  fitting 
the  curved  spectrum  of  a  source  with  a  power  law  (Abdo  et  al. 
2010a). 

Figure  16  shows  the  detection  efficiency  for  a  sample  like 
that  used  in  this  analysis.  Two  aspects  are  noteworthy.  First  the 
efficiency  at  Tfioo  =  10-8  photons  cm-2  s-1  (i.e.,  the  lowest 
flux  of  this  analysis  by  construction)  is  ~0.02  with  a  typical 
uncertainty  of  ±5  x  10-3  dictated  by  the  small  number  of 
sources  detected  in  our  simulations  at  the  lowest  fluxes.  Second, 
at  fluxes  around  Tqoo  ~  10-7  photons  cm-2  s-1  the  detection 
efficiency  becomes  larger  than  1.0.  This  effect  is  due  to  the  fact 


Figure  16.  Detection  efficiency  as  a  function  of  flux  for  a  population  of  sources 
with  curved  spectra  similar  to  those  of  FSRQs  determined  in  Section  5. 


that  fitting  a  curved  spectrum  source  with  a  power  law  yields  to 
an  overestimate  of  the  source  flux  by  a  factor  ~10%  (see  also 
Figure  8  in  Abdo  et  al.  201  Of).  Since  Figure  16  is  built  as  the 
ratio  (in  a  given  bin)  of  the  number  of  sources  detected  with  a 
given  flux  to  the  number  of  simulated  sources  with  that  flux,  the 
effect  mentioned  above  leads  to  a  detection  efficiency  >1.0. 

In  order  to  test  the  level  of  systematic  uncertainty  we  derived 
the  LF  using  the  detection  efficiency  reported  in  Figure  16 
(instead  of  using  the  detection  efficiency  curves  reported  in 
Abdo  et  al.  201  Of  and  used  throughout  Section  4).  Given  the 
“small”  number  of  sources  detected  in  the  three  simulations, 
it  was  not  possible  to  derive  a  two-dimensional  detection 
efficiency  as  a  function  of  flux  and  spectral  index  (like  that 
used  in  Section  4  and  derived  for  power-law  sources  in  Abdo 
et  al.  201  Of).  For  this  reason  the  parameters  of  the  distribution 
of  photon  indices  of  the  FSRQ  class  cannot  be  derived  from  the 
analysis  of  the  LF.  As  apparent  from  Table  3  most  parameters  of 
the  LF  derived  in  this  section  and  those  derived  in  Section  4.2 
are  compatible  within  their  statistical  errors.  The  only  parameter 
for  which  the  difference  is  slightly  larger  than  the  statistical 
errors  is  a.  The  parameter  a  governs  the  trend  of  the  redshift 
peak  with  luminosity  and  while  its  statistical  error  is  in  both 
case  0.03,  the  systematic  error  appears  to  be  0.05.  This  has 
very  little  impact  on  the  analysis  and  the  results  of  the  previous 
sections  are  fully  confirmed  and  robust  against  variations  of 
the  detection  efficiency  curve.  As  a  further  proof,  the  points  of 
the  de-evolved  LF  in  Figures  5  and  6  were  computed  using  the 
detection  efficiency  of  Figure  16  while  the  shaded  error  region 
was  computed  using  the  model  LF  derived  in  Section  4.2  that 
uses  the  detection  efficiency  for  power-law  sources. 

We  performed  an  additional  test  by  shifting  the  detection 
efficiency  curve  in  Figure  16  to  fluxes  10%  brighter  than 
measured.  The  rightward  shift  is  most  dramatic  as  it  increases  the 
magnitude  of  the  correction  at  faint  fluxes.  The  shift  is  performed 
in  order  to  account  for  uncertainties  in  the  determination  of  the 
detection  efficiency.  The  parameters  of  the  LF  are  all  consistent 
within  statistical  uncertainty  with  those  found  in  this  and  the 
previous  sections  and  reported  in  Table  3.  The  index  of  the  low- 
luminosity  slope  of  the  LF  becomes  slightly  steeper  (i.e.,  y\  = 
0.47  zb  0.18),  and  this  yields  a  slightly  larger  contribution  to 
the  IGRB  from  FSRQs.  We  thus  consider  the  typical  systematic 
uncertainty  connected  to  the  estimate  of  the  contribution  to  the 
IGRB  to  be  ~3%  of  the  IGRB  100  MeV-100  GeV  intensity. 
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A.  2.  Variability 

It  is  well  known  that  blazars  are  inherently  variable  objects 
with  variability  in  flux  of  up  to  a  factor  10  or  more.  Through¬ 
out  this  work  only  average  quantities  (i.e.,  mean  flux,  mean 
luminosity,  and  mean  photon  index)  are  used. 

It  is  not  straightforward  to  determine  how  blazar  variability 
affects  the  analysis  presented  here.  While  the  variability  patterns 
and  amplitudes  of  blazars  as  a  class  are  still  not  known  both 
Abdo  et  al.  (2010c)  and  Ackermann  et  al.  (2011)  presented 
a  detailed  analysis  of  the  variability  of  the  brightest  Fermi 
blazars.  They  report  that  the  variability  amplitude  of  the  FSRQ 
class  is  generally  larger  than  that  of  the  BL  Lac  population. 
However,  most  sources  (either  bright  or  faint  ones)  exceed  their 
average  flux  for  less  than  5%-20%  of  the  monitored  time  (i.e., 
respectively,  11  months  or  2  years).  This  drastically  reduces  the 
possibility  that  FSRQs  (or  blazars  more  in  general)  are  detected 
because  of  a  few  bright  flaring  episodes.  The  effect  of  high- 
amplitude  variability  connected  with  a  rising  density  of  sources 
at  smaller  flux  might  contaminate  samples,  as  the  one  used  here, 
with  objects  that  formally  should  not  be  included.  However, 
because  of  what  has  been  reported  above  and  the  flatness  of 
the  FSRQs  source  count  distribution  (see  Abdo  et  al.  201  Of; 
Ackermann  et  al.  2011,  and  Figure  2)  this  effect  is  very  likely 
marginal. 

Another,  smaller,  problem  is  connected  to  the  dependence  of 
the  effective  area  on  the  direction  of  the  incoming  photon  and  the 
LAT  detector  frame.16  Short  intense  flares  detected  during  fa¬ 
vorable  conditions  (i.e.,  on  axis  and  at  azimuthal  angles  of  ~0°, 
~90°,  ~180°,  or  ~270°)  might  lead  to  a  higher  TS,  increasing 
the  likelihood  of  source  detection.  However,  because  of  Fermi's 
continuous  scanning  of  the  sky  and  because  most  flares  are  ob¬ 
served  to  last  10  days  or  longer  (Abdo  et  al.  2010c),  the  effect 
above  has  negligible  influence  on  the  analysis  presented  here. 

Finally,  we  believe  that  variability  does  not  hamper  the  results 
of  this  analysis.  Using  even  longer  integration  times  (e.g.,  two 
or  three  years)  will  be  the  most  efficient  way  to  confirm  the 
results  of  this  analysis  and  dilute  the  effect  of  blazar  variability. 

A3.  Extragalactic  Background  Light 

Uncertainty  in  the  level  of  the  EBL,  in  particular  at  medium 
to  high  redshift,  might  in  principle  affect  our  analysis.  Energetic 
y- rays  from  FSRQs  at  high  redshifts  might  be  absorbed  by  the 
EBL  and  if  this  effect  is  not  taken  into  account  the  source-frame 
luminosity  would  be  underestimated.  This  would  lead  to  wrong 
estimates  of  the  space  densities  of  FSRQs.  However,  we  believe 
that  this  uncertainty  is  negligible. 

The  uncertainty  in  the  level  of  the  EBL  would  impact  the 
estimate  of  the  k-correction,  which  allows  us  to  determine  the 
source-frame  luminosities.  As  shown  in  Figure  10,  neglecting 
the  EBL  at  once  and  adopting  a  simple  power-law  spectrum  with 
a  photon  index  of  2.4  (instead  of  the  average  SED  determined  in 
Section  5)  introduces  an  uncertainty  of  ^4%  on  the  value  of  the 
k-correction  at  z  =  3.  Since  all  the  Fermi  FSRQs  are  detected 
within  this  redshift,  this  uncertainty  produces  no  appreciable 
impact  on  the  determination  of  the  LF. 
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